diff --git a/.gitignore b/.gitignore
index f1838d812a388a5b77ce0cb2f7cdbaa1ff3cee58..8e0f6676765f2e1de851784b336e74f176324381 100644
--- a/.gitignore
+++ b/.gitignore
@@ -46,4 +46,5 @@ out*
*.avi
gyacomo
gyacomo_dbg
-gyacomo23*
\ No newline at end of file
+gyacomo23*
+wk/parameters/profiles
\ No newline at end of file
diff --git a/matlab/compute/calculate_derivative.m b/matlab/compute/calculate_derivative.m
new file mode 100644
index 0000000000000000000000000000000000000000..5a37a900a3082723797c228eb9fc821b82219030
--- /dev/null
+++ b/matlab/compute/calculate_derivative.m
@@ -0,0 +1,30 @@
+function outputs = calculate_derivative(x, y, stencil_size)
+ % Check if the lengths of x and y are the same
+ if length(x) ~= length(y)
+ error('Input arrays x and y must have the same length.');
+ end
+
+ % Check if stencil_size is an odd number
+ if mod(stencil_size, 2) ~= 1
+ error('Stencil size must be an odd number for a central difference stencil.');
+ end
+
+ % Initialize the derivative array
+ derivative = zeros(size(y));
+
+ % Compute the derivative using central finite difference
+ half_stencil = (stencil_size - 1) / 2;
+ for i = (1 + half_stencil):(length(x) - half_stencil)
+ x_stencil = x(i - half_stencil : i + half_stencil);
+ y_stencil = y(i - half_stencil : i + half_stencil);
+
+ % Fit a polynomial and compute the derivative
+ p = polyfit(x_stencil, y_stencil, 1);
+ derivative(i) = p(1);
+ end
+
+ % Handle boundary points with forward/backward difference
+ derivative(1:half_stencil) = (y(half_stencil + 1) - y(1)) / (x(half_stencil + 1) - x(1));
+ derivative(end - half_stencil + 1:end) = (y(end) - y(end - half_stencil)) / (x(end) - x(end - half_stencil));
+ outputs = derivative;
+end
diff --git a/matlab/compute/compute_growth_rates.m b/matlab/compute/compute_growth_rates.m
index a71f8f061557bce683e54bef6783cb20e90a3db7..3d38c4bba2c0bbff9eec45b0bc75d0445e4efe89 100644
--- a/matlab/compute/compute_growth_rates.m
+++ b/matlab/compute/compute_growth_rates.m
@@ -1,4 +1,4 @@
-function [w,e,t] = compute_growth_rates(field,time)
+function [wavg, eavg, w, t] = compute_growth_rates(field,time)
% compute_growth_rates compute the linear growth rates using the amplitude
% ratio method
% see B.J. Frei et al. flux-tube paper
@@ -10,7 +10,8 @@ N1 = sz(1); N2 = sz(2); Nz = sz(3); Nt = sz(4);
w = zeros(N1,N2,Nt-1);
-e = zeros(N1,N2);
+wavg = zeros(N1,N2);
+eavg = zeros(N1,N2);
for i1 = 1:N1
for i2 = 1:N2
to_measure = reshape(field(i1,i2,:,:),Nz,Nt);
@@ -25,12 +26,17 @@ for i1 = 1:N1
w(i1,i2,it) = squeeze(sum(wl.*phi_nm1,1)./ZS);
end
% error estimation
- wavg = mean(w(i1,i2,ceil(Nt/2):end));
- wend = w(i1,i2,end);
- e(i1,i2) = abs(wend-wavg)/abs(wavg);
+ inputArray = real(squeeze(w(i1,i2,:)));
+ n = 5;
+ [avg_real, ~, sliceerr_real] = sliceAverage(inputArray, n);
+ inputArray = imag(squeeze(w(i1,i2,:)));
+ n = 5;
+ [avg_imag, ~, sliceerr_imag] = sliceAverage(inputArray, n);
+ wavg(i1,i2) = avg_real + 1i*avg_imag;
+ eavg(i1,i2) = mean(sliceerr_real) + 1i*mean(sliceerr_imag);
end
end
-t = time(2:Nt);
+t = time(1:Nt);
end
\ No newline at end of file
diff --git a/matlab/compute/compute_norm_grad.m b/matlab/compute/compute_norm_grad.m
new file mode 100644
index 0000000000000000000000000000000000000000..65efa80923fc1bcce1b235dd7caf0e5d36f92dd5
--- /dev/null
+++ b/matlab/compute/compute_norm_grad.m
@@ -0,0 +1,43 @@
+function [profile] = compute_norm_grad(profile_folder,xq)
+%COMPUTE_NORM_GRAD Compute the normalized gradients L_ref/L_chi=gradchi/chi
+ n_stencil = 3;
+
+ ne = load([profile_folder,'/ne.txt']);
+ x = ne(:,1);
+ y = ne(:,2);
+ y = interp1(x,y,xq);
+ % nmvm = ceil(numel(x)*smoothing_ratio);
+ % y = movmean(ne(:,2),nmvm);
+ % y = smooth_curve(x,ne(:,2),dx);
+
+ profile.Kne = calculate_derivative(xq,y,n_stencil);
+ profile.Kne(:,2) = profile.Kne(:,2)./y;
+ profile.ne = [xq y];
+
+ Te = load([profile_folder,'/Te.txt']);
+ x = Te(:,1);
+ y = Te(:,2);
+ y = interp1(x,y,xq);
+ % nmvm = ceil(numel(x)*smoothing_ratio);
+ % y = movmean(Te(:,2),nmvm);
+ % y = smooth_curve(x,Te(:,2),dx);
+ profile.KTe = calculate_derivative(xq,y,n_stencil);
+ profile.KTe(:,2) = profile.KTe(:,2)./y;
+ profile.Te = [xq y];
+
+ Ti = load([profile_folder,'/Ti.txt']);
+ x = Ti(:,1);
+ y = Ti(:,2);
+ y = interp1(x,y,xq);
+ % nmvm = ceil(numel(x)*smoothing_ratio);
+ % y = movmean(Ti(:,2),nmvm);
+ % y = smooth_curve(x,Ti(:,2),dx);
+ profile.KTi = calculate_derivative(xq,y,n_stencil);
+ profile.KTi(:,2) = profile.KTi(:,2)./y;
+ profile.Ti = [xq y];
+
+ % profile.nmvm= nmvm;
+ profile.n_stencil= n_stencil;
+ profile.xq = xq;
+ profile.fold=profile_folder;
+end
\ No newline at end of file
diff --git a/matlab/compute/interpolate_at_x0.m b/matlab/compute/interpolate_at_x0.m
new file mode 100644
index 0000000000000000000000000000000000000000..c7b170c0a20e1565749f58248c80909200ecd7a7
--- /dev/null
+++ b/matlab/compute/interpolate_at_x0.m
@@ -0,0 +1,9 @@
+function y0 = interpolate_at_x0(x, y, x0)
+ % Check if the lengths of x and y are the same
+ if length(x) ~= length(y)
+ error('Input arrays x and y must have the same length.');
+ end
+
+ % Perform linear interpolation using interp1
+ y0 = interp1(x, y, x0, 'linear', 'extrap');
+end
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index 0f04834b45e2c044c06720b014299c11ace6f31f..d449552ed0bc09b73e36b63f6241864e318e4abc 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -38,9 +38,13 @@ TW = [OPTIONS.KY_TW; OPTIONS.KX_TW];
[~,ikzf] = max(mean(squeeze(abs(field(1,:,FRAMES))),2));
% Amplitude ratio method to measure growth rates and freq.
-[wkykx, ekykx] = compute_growth_rates(FIELD(:,:,:,FRAMES),DATA.Ts3D(FRAMES));
+[~,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
FIGURE.fig = figure; %set(gcf, 'Position', [100 100 1200 700])
FIGURE.FIGNAME = 'mode_growth_meter';
for i = 1:2
@@ -54,17 +58,17 @@ for i = 1:2
case 'sum'
plt_ = @(f,ik) movmean(squeeze(sum(abs((f(:,ik,FRAMES))),1)),Nma);
MODESTR = '$\sum_{k_y}$';
- omega= @(ik) wkykx(1,:,end);
+ omega= @(ik) wkykx(1,:);
err = @(ik) ekykx(1,:);
case 'max'
plt_ = @(f,ik) movmean(squeeze(max(abs((f(:,ik,FRAMES))),[],1)),Nma);
MODESTR = '$\max_{k_y}$';
- omega= @(ik) wkykx(1,:,end);
+ omega= @(ik) wkykx(1,:);
err = @(ik) ekykx(1,:);
otherwise
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);
+ omega= @(ik) wkykx(OPTIONS.ik,:);
err = @(ik) ekykx(OPTIONS.ik,:);
end
kstr = 'k_x';
@@ -76,17 +80,17 @@ 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);
+ omega= @(ik) wkykx(:,1);
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);
+ omega= @(ik) wkykx(:,1);
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))];
- omega= @(ik) wkykx(:,OPTIONS.ik,end);
+ omega= @(ik) wkykx(:,OPTIONS.ik);
err = @(ik) ekykx(:,OPTIONS.ik);
end
kstr = 'k_y';
@@ -233,4 +237,5 @@ if d
end
top_title(DATA.paramshort)
+end
end
\ No newline at end of file
diff --git a/matlab/compute/smooth_curve.m b/matlab/compute/smooth_curve.m
new file mode 100644
index 0000000000000000000000000000000000000000..ea4280ad886213b25e89b72f479151810b0fdb08
--- /dev/null
+++ b/matlab/compute/smooth_curve.m
@@ -0,0 +1,15 @@
+function smoothed_y = smooth_curve(x, y, dx)
+ % Check if the lengths of x and y are the same
+ if length(x) ~= length(y)
+ error('Input arrays x and y must have the same length.');
+ end
+
+ % Initialize the smoothed_y array
+ smoothed_y = zeros(size(y));
+
+ % Perform smoothing
+ for i = 1:length(x)
+ indices = find(abs(x - x(i)) <= dx / 2);
+ smoothed_y(i) = mean(y(indices));
+ end
+end
\ No newline at end of file
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index 299c9b1119d5a8d2a7b5393db32b249fb05a2455..5d8f47a4c4d39ce199529ce938d989fb5a81bc0e 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -1,34 +1,35 @@
-function [FIG] = plot_ballooning(data,options)
+function [FIG, b_angle, phib, psib, ky_] = plot_ballooning(DATA,OPTIONS)
FIG.fig = figure; iplot = 1;
- [~,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));
+ [~,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));
ikyarray = unique(ikyarray);
- dz = data.grids.z(2) - data.grids.z(1);
- phi_real=real(data.PHI(:,:,:,it1));
- phi_imag=imag(data.PHI(:,:,:,it1));
+ dz = DATA.grids.z(2) - DATA.grids.z(1);
+ phi_real=real(DATA.PHI(:,:,:,it1));
+ phi_imag=imag(DATA.PHI(:,:,:,it1));
ncol = 1;
- if data.inputs.BETA > 0
- psi_real=real(data.PSI(:,:,:,it1));
- psi_imag=imag(data.PSI(:,:,:,it1));
+ if DATA.inputs.BETA > 0
+ psi_real=real(DATA.PSI(:,:,:,it1));
+ psi_imag=imag(DATA.PSI(:,:,:,it1));
ncol = 2;
end
- if options.PLOT_KP
+ if OPTIONS.PLOT_KP
ncol = 3;
end
% % Apply ballooning transform
- % if(data.grids.Nkx > 1)
- % nexc = round(data.grids.ky(2)*data.inputs.SHEAR*2*pi/data.grids.kx(2));
+ % if(DATA.grids.Nkx > 1)
+ % nexc = round(DATA.grids.ky(2)*DATA.inputs.SHEAR*2*pi/DATA.grids.kx(2));
% else
% nexc = 1;
% end
+ nline = 1;
for iky=ikyarray
- [phib_real,b_angle] = ballooning_representation(phi_real,data.inputs.SHEAR,data.grids.Npol,data.grids.kx,iky,data.grids.ky,data.grids.z);
- [phib_imag,~] = ballooning_representation(phi_imag,data.inputs.SHEAR,data.grids.Npol,data.grids.kx,iky,data.grids.ky,data.grids.z);
+ [phib_real,b_angle] = ballooning_representation(phi_real,DATA.inputs.SHEAR,DATA.grids.Npol,DATA.grids.kx,iky,DATA.grids.ky,DATA.grids.z);
+ [phib_imag,~] = ballooning_representation(phi_imag,DATA.inputs.SHEAR,DATA.grids.Npol,DATA.grids.kx,iky,DATA.grids.ky,DATA.grids.z);
phib = phib_real(:) + 1i * phib_imag(:);
% normalize real and imaginary parts at chi =0
- if options.normalized
+ if OPTIONS.normalized
[~,idxLFS] = min(abs(b_angle -0));
normalization = (phib( idxLFS));
else
@@ -39,24 +40,24 @@ function [FIG] = plot_ballooning(data,options)
phib_imag_norm = imag( phib_norm);
%
subplot(numel(ikyarray),ncol,ncol*(iplot-1)+1)
- plot(b_angle/pi, phib_real_norm,'-b'); hold on;
+ plot(b_angle/pi, phib_real_norm,'-b','DisplayName','$|\phi_B (\chi)|$'); hold on;
plot(b_angle/pi, phib_imag_norm ,'-r');
plot(b_angle/pi, sqrt(phib_real_norm .^2 + phib_imag_norm.^2),'-k');
legend('real','imag','norm')
xlabel('$\chi / \pi$')
ylabel('$\phi_B (\chi)$');
- title(['$k_y=',sprintf('%2.2f',data.grids.ky(iky)),...
- ',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']);
+ title(['$k_y=',sprintf('%2.2f',DATA.grids.ky(iky)),...
+ ',t_{avg}\in [',sprintf('%1.1f',DATA.Ts3D(it0)),',',sprintf('%1.1f',DATA.Ts3D(it1)),']$']);
% z domain start end point
- for i_ = 2*[1:data.grids.Nkx-1]-(data.grids.Nkx+1)
- xline(data.grids.Npol*(i_),'HandleVisibility','off'); hold on;
+ for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1)
+ xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on;
end
- for i_ = 2*[2:data.grids.Nkx]-(data.grids.Nkx+1)
- xline(data.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
+ for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1)
+ xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
end
- if data.inputs.BETA > 0
- [psib_real,b_angle] = ballooning_representation(psi_real,data.inputs.SHEAR,data.grids.Npol,data.grids.kx,iky,data.grids.ky,data.grids.z);
- [psib_imag,~] = ballooning_representation(psi_imag,data.inputs.SHEAR,data.grids.Npol,data.grids.kx,iky,data.grids.ky,data.grids.z);
+ if DATA.inputs.BETA > 0
+ [psib_real,b_angle] = ballooning_representation(psi_real,DATA.inputs.SHEAR,DATA.grids.Npol,DATA.grids.kx,iky,DATA.grids.ky,DATA.grids.z);
+ [psib_imag,~] = ballooning_representation(psi_imag,DATA.inputs.SHEAR,DATA.grids.Npol,DATA.grids.kx,iky,DATA.grids.ky,DATA.grids.z);
psib = psib_real(:) + 1i * psib_imag(:);
psib_norm = psib / normalization;
@@ -66,38 +67,50 @@ function [FIG] = plot_ballooning(data,options)
plot(b_angle/pi, psib_real_norm,'-b'); hold on;
plot(b_angle/pi, psib_imag_norm ,'-r');
plot(b_angle/pi, abs(psib_norm),'-k');
+
legend('real','imag','norm')
% z domain start end point
- for i_ = 2*[1:data.grids.Nkx-1]-(data.grids.Nkx+1)
- xline(data.grids.Npol*(i_),'HandleVisibility','off'); hold on;
+ for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1)
+ xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on;
end
- for i_ = 2*[2:data.grids.Nkx]-(data.grids.Nkx+1)
- xline(data.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
+ for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1)
+ xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
end
xlabel('$\chi / \pi$')
ylabel('$\psi_B (\chi)$');
- title(['$k_y=',sprintf('%2.2f',data.grids.ky(iky)),...
- ',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']);
+ title(['$k_y=',sprintf('%2.2f',DATA.grids.ky(iky)),...
+ ',t_{avg}\in [',sprintf('%1.1f',DATA.Ts3D(it0)),',',sprintf('%1.1f',DATA.Ts3D(it1)),']$']);
end
- if options.PLOT_KP
- kperp = h5read(data.outfilenames{1},'/data/grid/coordkp');
- [kperpb,b_angle] = ballooning_representation(kperp,data.inputs.SHEAR,data.grids.Npol,data.grids.kx,iky,data.grids.ky,data.grids.z);
+ if OPTIONS.PLOT_KP
+ kperp = h5read(DATA.outfilenames{1},'/DATA/grid/coordkp');
+ [kperpb,b_angle] = ballooning_representation(kperp,DATA.inputs.SHEAR,DATA.grids.Npol,DATA.grids.kx,iky,DATA.grids.ky,DATA.grids.z);
subplot(numel(ikyarray),ncol,ncol*(iplot-1)+3)
plot(b_angle/pi, kperpb,'-k'); hold on;
% z domain start end point
- for i_ = 2*[1:data.grids.Nkx-1]-(data.grids.Nkx+1)
- xline(data.grids.Npol*(i_),'HandleVisibility','off'); hold on;
+ for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1)
+ xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on;
end
- for i_ = 2*[2:data.grids.Nkx]-(data.grids.Nkx+1)
- xline(data.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
+ for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1)
+ xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
end
xlabel('$\chi / \pi$')
ylabel('$k_\perp (\chi)$');
- title(['$k_y=',sprintf('%2.2f',data.grids.ky(iky)),...
- ',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']);
+ title(['$k_y=',sprintf('%2.2f',DATA.grids.ky(iky)),...
+ ',t_{avg}\in [',sprintf('%1.1f',DATA.Ts3D(it0)),',',sprintf('%1.1f',DATA.Ts3D(it1)),']$']);
end
iplot = iplot + 1;
end
+ if ~OPTIONS.SHOWFIG
+ close
+ end
+ ky_ = DATA.grids.ky(iky);
end
+
+function [plots, ncurve] = add_to_plots(plots,plotname,x,xname,y,yname,ncurve)
+ plots.(plotname).x = x;
+ plots.(plotname).y = y;
+ plots.(plotname).xname = xname;
+ plots.(plotname).yname = yname;
+end
\ No newline at end of file
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
index 28e9b1c87cec3c969e6904bc1d3c52a4cf4e3d02..39d3855107422f1215465ced486610a4c5fc94a5 100644
--- a/matlab/plot/plot_metric.m
+++ b/matlab/plot/plot_metric.m
@@ -1,5 +1,5 @@
-function [ fig, arrays ] = 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'};
names = {'gxx','gxy','gxz','gyy','gyz','gzz',...
@@ -50,5 +50,13 @@ if NPLOT > 0
end
%outputs
arrays = squeeze(geo_arrays(:,:));
+R = [geo_arrays(:,12);geo_arrays(1,12)];
+Z = [geo_arrays(:,13);geo_arrays(1,13)];
+try
+ if ~options.SHOWFIG
+ close
+ end
+catch
+end
end
diff --git a/matlab/plot/plot_miller.m b/matlab/plot/plot_miller.m
new file mode 100644
index 0000000000000000000000000000000000000000..6577b064e2044f667871cb195be4df686ae64bb3
--- /dev/null
+++ b/matlab/plot/plot_miller.m
@@ -0,0 +1,18 @@
+function [R,Z] = plot_miller(geom,rho,Nz,PLOT)
+%
+R0 = geom.R0; Z0 = geom.Z0;
+kappa = geom.kappa; delta = geom.delta; zeta = geom.zeta;
+
+theta = linspace(0,2*pi,Nz+1);
+R = R0 + rho*cos(theta+asin(delta*sin(theta)));
+Z = Z0 + kappa*rho*sin(theta + zeta*sin(2*theta));
+
+if PLOT
+figure
+plot(R,Z,'-b');
+xlabel('R [m]');
+ylabel('Z [m]');
+axis tight
+axis equal
+end
+end
\ No newline at end of file
diff --git a/matlab/plot/show_geometry.m b/matlab/plot/show_geometry.m
index 01bc00d7c5bf341ce6f262b9d44bd380250eeb03..8cbd98eba6611ff6e4a4bd9bfdb87a8bf7d52fdc 100644
--- a/matlab/plot/show_geometry.m
+++ b/matlab/plot/show_geometry.m
@@ -3,7 +3,7 @@ function [ FIGURE ] = show_geometry(DATA,OPTIONS)
if DATA.grids.Nz > 1 % Torus flux tube geometry
Nturns = floor(abs(DATA.grids.z(1)/pi));
Nptor = ceil(DATA.grids.Nz*2/Nturns); Tgeom = 1;
- a = DATA.EPS; % Torus minor radius
+ a = DATA.inputs.EPS; % Torus minor radius
R = 1.; % Torus major radius
q = (DATA.grids.Nz>1)*DATA.inputs.Q0; % Flux tube safety factor
theta = linspace(-pi, pi, Nptor) ; % Poloidal angle
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 651998af5961619fb29f57c5a4769103b56ff72b..06290b4c8af281f769afe8aad7ff0da7f2b91bd9 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -238,7 +238,7 @@ switch OPTIONS.NAME
% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))),2));
tmp(:,:,iz) = abs((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))));
end
- FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
+ FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:)));
end
case 'Q_x' % ion heat flux
NAME = 'Qx';
@@ -254,7 +254,7 @@ switch OPTIONS.NAME
% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))),2));
tmp(:,:,iz) = abs((squeeze(fft2(qx_,DATA.Ny,DATA.Nx))));
end
- FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
+ FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:)));
end
case 'f_i'
SKIP_COMP = 1;
@@ -293,7 +293,7 @@ else
if REALP
tmp = zeros(Ny,Nx,Nz);
else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
+ tmp = zeros(DATA.grids.Nky,DATA.grids.Nkx,Nz);
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.grids.z)
diff --git a/matlab/setup.m b/matlab/setup.m
index 89d826810aec2bfb4792c01c2937a23d9799909b..7fd5e06d84ff5a1cfb74650a76a5e53afdf22e0b 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -82,8 +82,8 @@ switch CO
case 'LR'
COLL.mat_file = 'gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5';
case 'LD'
- % COLL.mat_file = 'gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5';
- COLL.mat_file = 'gk_landauii_P16_J9_dk_5e-2_km_2.0_NFLR_8.h5';
+ COLL.mat_file = 'gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5';
+ % COLL.mat_file = 'gk_landauii_P16_J9_dk_5e-2_km_2.0_NFLR_8.h5';
% COLL.mat_file = 'gk_landau_P11_J7_dk_5e-2_km_2.0_NFLR_16.h5';
% COLL.mat_file = 'gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
% COLL.mat_file = 'LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5';
diff --git a/wk/extract_miller_from_eqdsk.py b/wk/extract_miller_from_eqdsk.py
new file mode 100755
index 0000000000000000000000000000000000000000..e7138fd35c9bbc8857c77e0097d87d752b7cbd96
--- /dev/null
+++ b/wk/extract_miller_from_eqdsk.py
@@ -0,0 +1,645 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 26 09:32:25 2013
+
+@author: dtold
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from sys import exit, stdout
+import optparse as op
+
+nr = 120
+ntheta = 160
+parser = op.OptionParser(description='Extract Miller shaping parameters from EQDSK files.')
+parser.add_option('--rovera', '-r', action='store_const', const=1)
+parser.add_option('--conv', '-c', action='store_const', const=1)
+parser.add_option('--noplots', '-n', action='store_const', const=1)
+options, args = parser.parse_args()
+use_r_a = options.rovera
+show_plots = (not options.noplots)
+write_rhotor_rhopol_conversion = options.conv
+if not write_rhotor_rhopol_conversion:
+ if len(args) != 2:
+ exit("""
+Please give two arguments: <EQDSK filename> <Position in rho_tor>
+optional: -r <Position in r/a>
+ -c: Write rhotor/rhopol conversion to file
+ -n: Suppress plots\n""")
+
+filename = args[0]
+file = open(filename, 'r')
+
+
+def find(val, arr):
+ return np.argmin(np.abs(arr - val))
+
+
+eqdsk = file.readlines()
+# print('Header: {0:s}'.format(eqdsk[0]))
+# set resolutions
+nw = np.int(eqdsk[0].split()[-2])
+nh = np.int(eqdsk[0].split()[-1])
+pw = np.int((nw/8/2)*2) # psi-width, number of flux surfaces around position of interest
+# print('Resolution: {0:4d} x {1:4d}'.format(nw, nh))
+
+entrylength = 16
+try:
+ rdim, zdim, rctr, rmin, zmid = [np.float(eqdsk[1][j*entrylength:(j + 1)*entrylength]) for j in
+ range(len(eqdsk[1])//entrylength)]
+except:
+ entrylength = 15
+ try:
+ rdim, zdim, rctr, rmin, zmid = [np.float(eqdsk[1][j*entrylength:(j + 1)*entrylength]) for j
+ in range(len(eqdsk[1])//entrylength)]
+ except:
+ exit('Error reading EQDSK file, please check format!')
+rmag, zmag, psiax, psisep, Bctr = [np.float(eqdsk[2][j*entrylength:(j + 1)*entrylength]) for j in
+ range(len(eqdsk[2])//entrylength)]
+_, psiax2, _, rmag2, _ = [np.float(eqdsk[3][j*entrylength:(j + 1)*entrylength]) for j in
+ range(len(eqdsk[3])//entrylength)]
+zmag2, _, psisep2, _, _ = [np.float(eqdsk[4][j*entrylength:(j + 1)*entrylength]) for j in
+ range(len(eqdsk[4])//entrylength)]
+if rmag != rmag2:
+ np.sys.exit('Inconsistent rmag: %7.4g, %7.4g'%(rmag, rmag2))
+if psiax2 != psiax:
+ np.sys.exit('Inconsistent psiax: %7.4g, %7.4g'%(psiax, psiax2))
+if zmag != zmag2:
+ np.sys.exit('Inconsistent zmag: %7.4g, %7.4g'%(zmag, zmag2))
+if psisep2 != psisep:
+ np.sys.exit('Inconsistent psisep: %7.4g, %7.4g'%(psisep, psisep2))
+F = np.empty(nw, dtype=np.float)
+p = np.empty(nw, dtype=np.float)
+ffprime = np.empty(nw, dtype=np.float)
+pprime = np.empty(nw, dtype=np.float)
+qpsi = np.empty(nw, dtype=np.float)
+psirz_1d = np.empty(nw*nh, dtype=np.float)
+start_line = 5
+lines = range(nw//5)
+if nw%5 != 0:
+ lines = range(nw//5 + 1)
+for i in lines:
+ n_entries = len(eqdsk[i + start_line])//entrylength
+ F[i*5:i*5 + n_entries] = [np.float(eqdsk[i + start_line][j*entrylength:(j + 1)*entrylength]) for
+ j in range(n_entries)]
+start_line = i + start_line + 1
+
+for i in lines:
+ n_entries = len(eqdsk[i + start_line])//entrylength
+ p[i*5:i*5 + n_entries] = [np.float(eqdsk[i + start_line][j*entrylength:(j + 1)*entrylength]) for
+ j in range(n_entries)]
+start_line = i + start_line + 1
+
+for i in lines:
+ n_entries = len(eqdsk[i + start_line])//entrylength
+ ffprime[i*5:i*5 + n_entries] = [
+ np.float(eqdsk[i + start_line][j*entrylength:(j + 1)*entrylength]) for j in
+ range(n_entries)]
+start_line = i + start_line + 1
+
+for i in lines:
+ n_entries = len(eqdsk[i + start_line])//entrylength
+ pprime[i*5:i*5 + n_entries] = [
+ np.float(eqdsk[i + start_line][j*entrylength:(j + 1)*entrylength]) for j in
+ range(n_entries)]
+start_line = i + start_line + 1
+
+lines_twod = range(nw*nh//5)
+if nw*nh%5 != 0:
+ lines_twod = range(nw*nh//5 + 1)
+for i in lines_twod:
+ n_entries = len(eqdsk[i + start_line])//entrylength
+ psirz_1d[i*5:i*5 + n_entries] = [
+ np.float(eqdsk[i + start_line][j*entrylength:(j + 1)*entrylength]) for j in
+ range(n_entries)]
+start_line = i + start_line + 1
+psirz = psirz_1d.reshape(nh, nw)
+
+for i in lines:
+ n_entries = len(eqdsk[i + start_line])//entrylength
+ qpsi[i*5:i*5 + n_entries] = [np.float(eqdsk[i + start_line][j*entrylength:(j + 1)*entrylength])
+ for j in range(n_entries)]
+start_line = i + start_line + 1
+
+# invert sign of psi if necessary to guarantee increasing values for interpolation
+if psisep < psiax:
+ psirz = -psirz
+ ffprime = -ffprime
+ pprime = -pprime
+ psiax *= -1
+ psisep *= -1
+
+# ignore limiter data etc. for the moment
+dw = rdim/(nw - 1)
+dh = zdim/(nh - 1)
+rgrid = np.array([rmin + i*dw for i in range(nw)])
+zgrid = np.array([zmid - zdim/2. + i*dh for i in range(nh)])
+# contourf(rgrid,zgrid,psirz,70);gca().set_aspect('equal')
+# show()
+
+# create 5th order 2D spline representation of Psi(R,Z)
+from scipy.interpolate import RectBivariateSpline
+from scipy.interpolate import interp1d
+from scipy.interpolate import UnivariateSpline
+
+interpol_order = 3
+psi_spl = RectBivariateSpline(zgrid, rgrid, psirz, kx=interpol_order, ky=interpol_order)
+
+# linear grid of psi, on which all 1D fields are defined
+linpsi = np.linspace(psiax, psisep, nw)
+# create rho_tor grid
+x_fine = np.linspace(psiax, psisep, nw*10)
+phi_fine = np.empty((nw*10), dtype=np.float)
+
+phi_fine[0] = 0.
+# spline of q for rhotor grid
+q_spl_psi = UnivariateSpline(linpsi, qpsi, k=interpol_order, s=1e-5)
+q_fine = q_spl_psi(x_fine)
+
+for i in range(1, nw*10):
+ x = x_fine[:i + 1]
+ y = q_spl_psi(x)
+ phi_fine[i] = np.trapz(y, x)
+rho_tor_fine = np.sqrt(phi_fine/phi_fine[-1])
+rho_tor_spl = UnivariateSpline(x_fine, rho_tor_fine, k=interpol_order, s=1e-5)
+rho_tor = np.empty(nw, dtype=np.float)
+for i in range(nw):
+ rho_tor[i] = rho_tor_spl(linpsi[i])
+
+if write_rhotor_rhopol_conversion:
+ rt_rp_filename = 'rt_rp_%s'%filename
+ rt_rp_file = open(rt_rp_filename, 'w')
+ rt_rp_file.write('# rho_tor rho_pol q\n')
+ for i in range(len(x_fine)):
+ rho_pol = np.sqrt((x_fine[i] - psiax)/(psisep - psiax))
+ rt_rp_file.write('%16.8e %16.8e %16.8e\n'%(rho_tor_fine[i], rho_pol, q_fine[i]))
+ rt_rp_file.close()
+ exit('\nWrote rhotor/rhopol conversion to %s.'%rt_rp_filename)
+
+t1 = np.arctan2(zmid - zdim/2. - zmag, rmin - rmag)
+t2 = np.arctan2(zmid - zdim/2. - zmag, rmin + rdim - rmag)
+t3 = np.arctan2(zmid + zdim/2. - zmag, rmin + rdim - rmag)
+t4 = np.arctan2(zmid + zdim/2. - zmag, rmin - rmag)
+
+theta_arr = np.linspace(-np.pi, np.pi, ntheta)
+# for i in range(nw):
+# curr_psi=linpsi[i]
+
+# print('Finding flux surface shapes...')
+R = np.empty((nw, ntheta), dtype=np.float)
+Z = np.empty((nw, ntheta), dtype=np.float)
+dr = rdim*np.cos(theta_arr)
+dz = rdim*np.sin(theta_arr)
+
+for j in range(len(theta_arr)):
+ # stdout.write('\r Finished %4.1f%%.'%(j*100./(ntheta - 1)))
+ stdout.flush()
+ theta = theta_arr[j]
+ r_pol = np.linspace(rmag, rmag + dr[j], nr) # array([rmag+i*dr for i in range(nr)])
+ z_pol = np.linspace(zmag, zmag + dz[j], nr) # array([zmag+i*dz for i in range(nr)])
+ psi_rad = psi_spl.ev(z_pol, r_pol)
+ psi_rad_sav = psi_rad
+ psi_rad[0] = psiax
+ # must restrict interpolation range because of non-monotonic psi around coils
+ cutoff = 0
+ for i in range(1, len(psi_rad)):
+ if psi_rad[i] < psi_rad[i - 1]:
+ cutoff = i
+ break
+ psi_rad = psi_rad[:i]
+ end_ind = np.argmin(np.abs(psi_rad - psisep))
+ end_ind += (1 if (psi_rad[end_ind] < psisep) else 0)
+ indsep = end_ind + 1
+
+ R_int = interp1d(psi_rad[:indsep], r_pol[:indsep], kind=interpol_order)
+ R[:, j] = R_int(linpsi)
+ Z_int = interp1d(psi_rad[:indsep], z_pol[:indsep], kind=interpol_order)
+ Z[:, j] = Z_int(linpsi)
+
+# print('\nFinding flux surface centers...')
+# find average elevation for all flux surfaces
+Z_avg = np.empty(nw, dtype=np.float)
+ds = np.empty(ntheta, dtype=np.float)
+for i in range(1, nw):
+ ds[1:ntheta - 1] = 0.5*np.sqrt(
+ (R[i, 2:ntheta] - R[i, 0:ntheta - 2]) ** 2 + (Z[i, 2:ntheta] - Z[i, 0:ntheta - 2]) ** 2)
+ ds[0] = 0.5*np.sqrt((R[i, 1] - R[i, -1]) ** 2 + (Z[i, 1] - Z[i, -1]) ** 2)
+ ds[-1] = 0.5*np.sqrt((R[i, 0] - R[i, -2]) ** 2 + (Z[i, 0] - Z[i, -2]) ** 2)
+ Z_avg[i] = np.average(Z[i, :], weights=ds)
+# Treat the magnetic axis separately as no average is required and ds==0
+Z_avg[0] = Z[0, 0]
+
+# find R0,Z0 for all flux surfaces
+R0 = np.empty(nw, dtype=np.float)
+R0[0] = rmag
+r_avg = np.empty(nw, dtype=np.float)
+r_avg[0] = 0.
+r_maxmin = np.empty(nw, dtype=np.float)
+r_maxmin[0] = 0.
+for i in range(1, nw):
+ # stdout.write('\r Finished %4.1f%%.'%(i*100./(nw - 1)))
+ stdout.flush()
+ R_array = R[i, ntheta//4:3*ntheta//4]
+ Z_array = Z[i, ntheta//4:3*ntheta//4]
+ # low field side
+ Z_int = interp1d(Z_array, range(ntheta//2), kind=interpol_order)
+ ind_Zavg = Z_int(Z_avg[i])
+ R_int = interp1d(range(ntheta//2), R_array, kind=interpol_order)
+ R_out = R_int(ind_Zavg)
+ R_max = np.amax(R_array)
+ # high field side
+ R_array = np.roll(R[i, :-1], ntheta//2)[ntheta//4:3*ntheta//4]
+ Z_array = np.roll(Z[i, :-1], ntheta//2)[ntheta//4:3*ntheta//4]
+
+ # have to use negative Z_array here to have increasing order
+ Z_int = interp1d(-Z_array, range(ntheta//2), kind=interpol_order)
+ # again negative
+ ind_Zavg = Z_int(-Z_avg[i])
+
+ R_int = interp1d(range(ntheta//2), R_array, kind=interpol_order)
+ R_in = R_int(ind_Zavg)
+ R_min = np.amin(R_array)
+ R0[i] = 0.5*(R_out + R_in)
+ r_avg[i] = 0.5*(R_out - R_in)
+ r_maxmin[i] = 0.5*(R_max - R_min)
+
+radpos = np.float(args[1])
+if use_r_a:
+ r_a = radpos
+ # find psi index of interest (for the specified r/a position)
+ poi_ind = find(radpos, r_avg/r_avg[-1])
+ # print('\nExamine {0:3d} flux surfaces around position rho_tor={1:7.4g}...'.format(pw, radpos))
+else:
+ # find psi index of interest (for the specified rho_tor position)
+ poi_ind = find(radpos, rho_tor)
+ ravg_rho_spl = UnivariateSpline(rho_tor, r_avg/r_avg[-1], k=interpol_order, s=1e-5)
+ r_a = np.float(ravg_rho_spl(radpos))
+ # print('\nExamine {0:3d} flux surfaces around position r/a={1:7.4g}...'.format(pw, r_a))
+
+# modified theta grid for each flux surface
+# arrays equidistant on modified theta grid are marked by 'tm' index!!!
+linpsi_spl = UnivariateSpline(r_avg, linpsi, k=interpol_order, s=1e-5)
+ravg_spl = UnivariateSpline(linpsi, r_avg, k=interpol_order, s=1e-5)
+
+rmaxmin_spl = UnivariateSpline(linpsi, r_maxmin, k=interpol_order, s=1e-5)
+q_spl = UnivariateSpline(r_avg, qpsi, k=interpol_order, s=1e-5)
+R0_spl = UnivariateSpline(r_avg, R0, k=interpol_order, s=1e-5)
+Z0_spl = UnivariateSpline(r_avg, Z_avg, k=interpol_order, s=1e-5)
+F_spl = UnivariateSpline(r_avg, F, k=interpol_order, s=1e-5)
+r = r_a*r_avg[-1]
+psi = np.float(linpsi_spl(r))
+psi_N = (psi - psiax)/(psisep - psiax)
+R0_pos = np.float(R0_spl(r))
+Z0_pos = np.float(Z0_spl(r))
+F_pos = np.float(F_spl(r))
+Bref_miller = F_pos/R0_pos
+
+# print('Coordinates: r={0:8.5g}, psi={1:8.5g}, psi_N={2:8.5g}, r/R0={3:8.5g}, rho_tor={4:8.5g}, '
+ # 'r_maxmin={5:8.5g}'.format(r, psi, psi_N, r/R0_pos, np.float(rho_tor_spl(psi)),
+ # np.float(rmaxmin_spl(psi))))
+psi_stencil = range(poi_ind - pw//2, poi_ind + pw//2)
+if psi_stencil[0] < 1:
+ psi_stencil = [psi_stencil[i] + 1 - psi_stencil[0] for i in range(len(psi_stencil))]
+if psi_stencil[-1] > nw - 1:
+ psi_stencil = [psi_stencil[i] - (psi_stencil[-1] - nw + 1) for i in range(len(psi_stencil))]
+R_tm = np.empty((pw, ntheta), dtype=np.float)
+Z_tm = np.empty((pw, ntheta), dtype=np.float)
+R_extended = np.empty(2*ntheta - 1, dtype=np.float)
+Z_extended = np.empty(2*ntheta - 1, dtype=np.float)
+# theta_mod[0]=theta_arr
+# R_tm[0]=R[0]
+# Z_tm[0]=Z[0]
+theta_tmp = np.linspace(-2.*np.pi, 2*np.pi, 2*ntheta - 1)
+
+# print('Interpolating to flux-surface dependent (proper) theta grid...')
+for i in psi_stencil:
+ # stdout.write('\r Finished %4.1f%%.'%(psi_stencil.index(i)*100./(len(psi_stencil) - 1)))
+ stdout.flush()
+ imod = i - psi_stencil[0]
+ # print('Finished {0:4.1f}%%.'.format(float(i)/(pw-1)*100))
+ R_extended[0:(ntheta - 1)//2] = R[i, (ntheta + 1)//2:-1]
+ R_extended[(ntheta - 1)//2:(3*ntheta - 3)//2] = R[i, :-1]
+ R_extended[(3*ntheta - 3)//2:] = R[i, 0:(ntheta + 3)//2]
+ Z_extended[0:(ntheta - 1)//2] = Z[i, (ntheta + 1)//2:-1]
+ Z_extended[(ntheta - 1)//2:(3*ntheta - 3)//2] = Z[i, :-1]
+ Z_extended[(3*ntheta - 3)//2:] = Z[i, 0:(ntheta + 3)//2]
+ # for j in range(ntheta):
+ theta_mod_ext = np.arctan2(Z_extended - Z_avg[i], R_extended - R0[i])
+ # introduce 2pi shifts to theta_mod_ext
+ for ind in range(ntheta):
+ if theta_mod_ext[ind + 1] < 0. and theta_mod_ext[ind] > 0. and np.abs(
+ theta_mod_ext[ind + 1] - theta_mod_ext[ind]) > np.pi:
+ lshift_ind = ind
+ if theta_mod_ext[-ind - 1] > 0. and theta_mod_ext[-ind] < 0. and np.abs(
+ theta_mod_ext[-ind - 1] - theta_mod_ext[-ind]) > np.pi:
+ rshift_ind = ind
+ theta_mod_ext[-rshift_ind:] += 2.*np.pi
+ theta_mod_ext[:lshift_ind + 1] -= 2.*np.pi
+ # print theta_mod, theta_arr
+ # plot(theta_mod_ext)
+ # plot(theta_tmp)
+ # show()
+ theta_int = interp1d(theta_mod_ext, theta_tmp, kind=interpol_order)
+ theta_orig_tm = theta_int(theta_arr)
+ R_int = interp1d(theta_mod_ext, R_extended, kind=interpol_order)
+ Z_int = interp1d(theta_mod_ext, Z_extended, kind=interpol_order)
+ R_tm[imod] = R_int(theta_arr)
+ Z_tm[imod] = Z_int(theta_arr)
+# plot(R_tm[imod],Z_tm[imod])
+# gca().set_aspect('equal')
+
+# now we have the flux surfaces on a symmetric grid in theta (with reference to R0(r), Z0(r))
+# symmetrize flux surfaces
+# figure()
+R_sym = np.empty((pw, ntheta), dtype=np.float)
+Z_sym = np.empty((pw, ntheta), dtype=np.float)
+for i in psi_stencil:
+ imod = i - psi_stencil[0]
+ Z_sym[imod, :] = 0.5*(Z_tm[imod, :] - Z_tm[imod, ::-1]) + Z_avg[i]
+ R_sym[imod, :] = 0.5*(R_tm[imod, :] + R_tm[imod, ::-1])
+# plot(R_sym[imod],Z_sym[imod])
+# gca().set_aspect('equal')
+# show()
+
+dq_dr_avg = np.empty(pw, dtype=np.float)
+dq_dpsi = np.empty(pw, dtype=np.float)
+drR = np.empty(pw, dtype=np.float)
+drZ = np.empty(pw, dtype=np.float)
+kappa = np.empty(pw, dtype=np.float)
+delta = np.empty(pw, dtype=np.float)
+s_kappa = np.empty(pw, dtype=np.float)
+s_delta = np.empty(pw, dtype=np.float)
+delta_upper = np.empty(pw, dtype=np.float)
+delta_lower = np.empty(pw, dtype=np.float)
+zeta_arr = np.empty((pw, 4), dtype=np.float)
+zeta = np.empty(pw, dtype=np.float)
+s_zeta = np.empty(pw, dtype=np.float)
+for i in psi_stencil:
+ imod = i - psi_stencil[0]
+ # calculate delta
+ stencil_width = ntheta//10
+ for o in range(2):
+ if o:
+ ind = np.argmax(Z_sym[imod])
+ section = range(ind + stencil_width//2, ind - stencil_width//2, -1)
+ else:
+ ind = np.argmin(Z_sym[imod])
+ section = range(ind - stencil_width//2, ind + stencil_width//2)
+ x = R_sym[imod, section]
+ y = Z_sym[imod, section]
+ y_int = interp1d(x, y, kind=interpol_order)
+ x_fine = np.linspace(np.amin(x), np.amax(x), stencil_width*100)
+ y_fine = y_int(x_fine)
+ if o:
+ x_at_extremum = x_fine[np.argmax(y_fine)]
+ delta_upper[imod] = (R0[i] - x_at_extremum)/r_avg[i]
+ Z_max = np.amax(y_fine)
+ else:
+ x_at_extremum = x_fine[np.argmin(y_fine)]
+ delta_lower[imod] = (R0[i] - x_at_extremum)/r_avg[i]
+ Z_min = np.amin(y_fine)
+ # calculate kappa
+ kappa[imod] = (Z_max - Z_min)/2./r_avg[i]
+
+# linear extrapolation (in psi) for axis values
+# delta_upper[0]=2*delta_upper[1]-delta_upper[2]
+# delta_lower[0]=2*delta_lower[1]-delta_lower[2]
+# kappa[0]=2*kappa[1]-kappa[2]
+# zeta[0]=2*zeta[1]-zeta[2]
+delta = 0.5*(delta_upper + delta_lower)
+
+# calculate zeta
+for i in psi_stencil:
+ imod = i - psi_stencil[0]
+ x = np.arcsin(delta[imod])
+ # find the points that correspond to Miller-theta=+-pi/4,+-3/4*pi and extract zeta from those
+ for o in range(4):
+ if o == 0:
+ val = np.pi/4.
+ searchval = np.cos(val + x/np.sqrt(2))
+ searcharr = (R_sym[imod] - R0[i])/r_avg[i]
+ elif o == 1:
+ val = 3.*np.pi/4
+ searchval = np.cos(val + x/np.sqrt(2))
+ searcharr = (R_sym[imod] - R0[i])/r_avg[i]
+ elif o == 2:
+ val = -np.pi/4.
+ searchval = np.cos(val - x/np.sqrt(2))
+ searcharr = (R_sym[imod] - R0[i])/r_avg[i]
+ elif o == 3:
+ val = -3.*np.pi/4
+ searchval = np.cos(val - x/np.sqrt(2))
+ searcharr = (R_sym[imod] - R0[i])/r_avg[i]
+ if o in [0, 1]:
+ searcharr2 = searcharr[ntheta//2:]
+ ind = find(searchval, searcharr2) + ntheta//2
+ else:
+ searcharr2 = searcharr[0:ntheta//2]
+ ind = find(searchval, searcharr2)
+ # print o,ind
+ section = range(ind - stencil_width//2, ind + stencil_width//2)
+ theta_sec = theta_arr[section]
+ if o in [0, 1]:
+ theta_int = interp1d(-searcharr[section], theta_sec, kind=interpol_order)
+ theta_of_interest = theta_int(-searchval)
+ else:
+ theta_int = interp1d(searcharr[section], theta_sec, kind=interpol_order)
+ theta_of_interest = theta_int(searchval)
+ Z_sec = Z_sym[imod, section]
+ Z_sec_int = interp1d(theta_sec, Z_sec, kind=interpol_order)
+ # print searchval,val, theta_sec
+ Z_val = Z_sec_int(theta_of_interest)
+ zeta_arr[imod, o] = np.arcsin((Z_val - Z_avg[i])/kappa[imod]/r_avg[i])
+ zeta_arr[imod, 1] = np.pi - zeta_arr[imod, 1]
+ zeta_arr[imod, 3] = -np.pi - zeta_arr[imod, 3]
+ # print zeta_arr[i]
+ zeta[imod] = 0.25*(
+ np.pi + zeta_arr[imod, 0] - zeta_arr[imod, 1] - zeta_arr[imod, 2] + zeta_arr[imod, 3])
+
+Bref_efit = np.abs(F[0]/R0[0])
+Lref_efit = np.sqrt(2*np.abs(phi_fine[-1])/Bref_efit)
+
+kappa_spl = UnivariateSpline(r_avg[psi_stencil], kappa, k=interpol_order, s=1e-5)
+delta_spl = UnivariateSpline(r_avg[psi_stencil], delta, k=interpol_order, s=1e-5)
+zeta_spl = UnivariateSpline(r_avg[psi_stencil], zeta, k=interpol_order, s=1e-5)
+amhd = np.empty(pw, dtype=np.float)
+amhd_Miller = np.empty(pw, dtype=np.float)
+Vprime = np.empty(pw, dtype=np.float)
+dV_dr = np.empty(pw, dtype=np.float)
+V = np.empty(pw, dtype=np.float)
+V_manual = np.empty(pw, dtype=np.float)
+r_FS = np.empty(pw, dtype=np.float)
+for i in psi_stencil:
+ imod = i - psi_stencil[0]
+ Vprime[imod] = np.abs(np.sum(qpsi[i]*R_sym[imod] ** 2/F[i])*4*np.pi ** 2/ntheta)
+ dV_dr[imod] = np.abs(np.sum(qpsi[i]*R_sym[imod] ** 2/F[i])*4*np.pi ** 2/ntheta)/ \
+ ravg_spl.derivatives(linpsi[i])[1]
+ # V[imod]=trapz(Vprime[:imod+1],linpsi[psi_stencil])
+ r_FS[imod] = np.average(np.sqrt((R_sym[imod] - R0[i]) ** 2 + (Z_sym[imod] - Z_avg[i]) ** 2),
+ weights=qpsi[i]*R_sym[imod] ** 2/F[i])
+ amhd[imod] = -qpsi[i] ** 2*R0[i]*pprime[i]*8*np.pi*1e-7/Bref_miller ** 2/ \
+ ravg_spl.derivatives(linpsi[i])[1]
+ # amhd_Miller[imod]=-2*Vprime[imod]/(2*pi)**2*(V[imod]/2/pi**2/R0[i])**0.5*4e-7*pi*pprime[i]
+ dq_dr_avg[imod] = q_spl.derivatives(r_avg[i])[1]
+ dq_dpsi[imod] = q_spl_psi.derivatives(linpsi[i])[1]
+ drR[imod] = R0_spl.derivatives(r_avg[i])[1]
+ drZ[imod] = Z0_spl.derivatives(r_avg[i])[1]
+ s_kappa[imod] = kappa_spl.derivatives(r_avg[i])[1]*r_avg[i]/kappa[imod]
+ s_delta[imod] = delta_spl.derivatives(r_avg[i])[1]*r_avg[i]/np.sqrt(1 - delta[imod] ** 2)
+ s_zeta[imod] = zeta_spl.derivatives(r_avg[i])[1]*r_avg[i]
+amhd_spl = UnivariateSpline(r_avg[psi_stencil], amhd, k=interpol_order, s=1e-5)
+rFS_spl = UnivariateSpline(r_avg[psi_stencil], r_FS, k=interpol_order, s=1e-5)
+drR_spl = UnivariateSpline(r_avg[psi_stencil], drR, k=interpol_order, s=1e-5)
+drZ_spl = UnivariateSpline(r_avg[psi_stencil], drZ, k=interpol_order, s=1e-5)
+Zavg_spl = UnivariateSpline(r_avg, Z_avg, k=interpol_order, s=1e-5)
+
+# plot(r_avg[psi_stencil],V)
+# figure()
+fig = plt.figure()
+plt.rcParams.update({'axes.titlesize': 'small',
+ 'axes.labelsize': 'small',
+ 'xtick.labelsize': 'small',
+ 'ytick.labelsize': 'small',
+ 'legend.fontsize': 'small'})
+
+ax = fig.add_subplot(3, 3, 2)
+ax.plot(r_avg[psi_stencil], kappa)
+ax.set_title('Elongation')
+ax.set_xlabel(r'$r_{avg}$')
+ax.set_ylabel(r'$\kappa$')
+ax.axvline(r, 0, 1, ls='--', color='k', lw=2)
+
+ax2 = fig.add_subplot(3, 3, 3)
+ax2.plot(r_avg[psi_stencil], s_kappa)
+ax2.set_title('Elongation (Derivative)')
+ax2.set_xlabel(r'$r_{avg}$')
+ax2.set_ylabel(r'$s_\kappa$')
+ax2.axvline(r, 0, 1, ls='--', color='k', lw=2)
+
+ax3 = fig.add_subplot(3, 3, 5)
+ax3.plot(r_avg[psi_stencil], delta)
+ax3.set_title('Triangularity')
+ax3.set_xlabel(r'$r_{avg}$')
+ax3.set_ylabel(r'$\delta$')
+ax3.axvline(r, 0, 1, ls='--', color='k', lw=2)
+
+ax4 = fig.add_subplot(3, 3, 6)
+ax4.plot(r_avg[psi_stencil], s_delta)
+ax4.set_title('Triangularity (Derivative)')
+ax4.set_xlabel(r'$r_{avg}$')
+ax4.set_ylabel(r'$s_\delta$')
+ax4.axvline(r, 0, 1, ls='--', color='k', lw=2)
+# figure()
+# plot(r_avg[psi_stencil],delta_diff)
+# xlabel(r'$r_{avg}$',fontsize=14)
+# ylabel(r'$\delta_u-\delta_l$',fontsize=14)
+
+ax5 = fig.add_subplot(3, 3, 8)
+ax5.plot(r_avg[psi_stencil], zeta)
+ax5.set_title('Squareness')
+ax5.set_xlabel(r'$r_{avg}$')
+ax5.set_ylabel(r'$\zeta$')
+ax5.axvline(r, 0, 1, ls='--', color='k', lw=2)
+
+ax6 = fig.add_subplot(3, 3, 9)
+ax6.plot(r_avg[psi_stencil], s_zeta)
+ax6.set_title('Squareness (Derivative)')
+ax6.set_xlabel(r'$r_{avg}$')
+ax6.set_ylabel(r'$s_\zeta$')
+ax6.axvline(r, 0, 1, ls='--', color='k', lw=2)
+
+# select a given flux surface
+ind = poi_ind - psi_stencil[0] # find(r_a,r_avg/r_avg[-1])
+Lref = R0_pos
+# print('\n\nShaping parameters for flux surface r={0:9.5g}, r/a={1:9.5g}:'.format(r, r_a))
+# print 'r_FS= %9.5g (flux-surface averaged radius)\n' %rFS_spl(r)
+# print('Copy the following block into a GENE parameters file:\n')
+print('trpeps = {0:9.5g}'.format(r/R0_pos))
+print('q0 = {0:9.5g}'.format(np.float(q_spl(r))))
+print('shat = {0:9.5g} '.format(
+ r/np.float(q_spl(r))*q_spl.derivatives(r)[1]))
+# print 'shat=%9.5g (defined as (psi-psiax)/q*dq_dpsi)' %((psi-linpsi[0])/q_spl(
+# r)*q_spl_psi.derivatives(psi)[1])
+print('amhd = {0:9.5g}'.format(np.float(amhd_spl(r))))
+# print 'amhd_Miller=%9.5g' %amhd_Miller[ind]
+# print test[ind]
+print('drR = {0:9.5g}'.format(np.float(drR_spl(r))))
+print('drZ = {0:9.5g}'.format(np.float(drZ_spl(r))))
+print('kappa = {0:9.5g}'.format(np.float(kappa_spl(r))))
+print('s_kappa = {0:9.5g}'.format(kappa_spl.derivatives(r)[1]*r/np.float(kappa_spl(r))))
+print('delta = {0:9.5g}'.format(np.float(delta_spl(r))))
+print('s_delta = {0:9.5g}'.format(
+ delta_spl.derivatives(r)[1]*r/np.sqrt(1 - np.float(delta_spl(r)) ** 2)))
+print('zeta = {0:9.5g}'.format(np.float(zeta_spl(r))))
+print('s_zeta = {0:9.5g}'.format(zeta_spl.derivatives(r)[1]*r))
+print('minor_r = {0:9.5g}'.format(1.0))
+print('major_R = {0:9.5g}'.format(R0_pos/r*r_a))
+# print('\nFor normalization to major radius, set instead:')
+print('minor_r = {0:9.5g}'.format(r/R0_pos/r_a))
+print('major_R = {0:9.5g}'.format(1.0))
+# print('The same conversion factor must be considered for input frequencies and gradients.')
+# print('\nAdditional information:')
+# print('Lref = {0:9.5g} !for Lref=a convention'.format(r_avg[-1]))
+print('Lref = {0:9.5g}'.format(R0_pos))
+print('Bref = {0:9.5g}'.format(Bref_miller))
+# minor radius at average flux surface elevation (GYRO definition)
+print('a = {0:9.5g}'.format(r_avg[-1]))
+print('R0 = {0:9.5g}'.format(R0_pos))
+print('Z0 = {0:9.5g}'.format(np.float(Zavg_spl(r))))
+print('Lref_efit = {0:9.5g}'.format(Lref_efit))
+print('Bref_efit = {0:9.5g}'.format(Bref_efit))
+print('B_unit_GYRO= {0:9.5g}'.format(q_spl(r)/r/ravg_spl.derivatives(psi)[1]))
+# print 'Vprime=%9.5g; dV_dr=%9.5g' %(Vprime[ind],dV_dr[ind])
+# print 'V=%9.5g' %V[ind]
+# minor radius defined by 0.5(R_max-R_min), where R_max and R_min can have any elevation
+# print 'a_maxmin = %9.5g' %r_maxmin[-1]
+print('dpsidr = {0:9.5g}'.format(1./ravg_spl.derivatives(psi)[1]))
+# print 'dr_maxmin/dr = %9.5g' %(rmaxmin_spl.derivatives(psi)[1]/ravg_spl.derivatives(psi)[1])
+print('drho_tordr = {0:9.5g}'.format(rho_tor_spl.derivatives(psi)[1]/ravg_spl.derivatives(psi)[1]))
+# print('Gradient conversion omt(rho_tor) -> a/LT; factor = {0:9.5g}'.format(
+ # r_avg[-1]*(rho_tor_spl.derivatives(psi)[1]/ravg_spl.derivatives(psi)[1])))
+# print('Gradient conversion omt(rho_tor) -> R/LT; factor = {0:9.5g}'.format(
+ # R0_pos*(rho_tor_spl.derivatives(psi)[1]/ravg_spl.derivatives(psi)[1])))
+
+#fig7 = plt.figure()
+ax7 = fig.add_subplot(1, 3, 1)
+ax7.plot(R_tm[ind], Z_tm[ind], 'k-', lw=2, label='original')
+ax7.plot(R_sym[ind], Z_sym[ind], 'r-', lw=2, label='symmetrized')
+# for v in range(-10,10,4):
+# zeta_test=-0.011086#zeta[ind]
+ax7.plot(R0_pos + r*np.cos(theta_arr + np.arcsin(delta_spl(r))*np.sin(theta_arr)),
+ Zavg_spl(r) + kappa_spl(r)*r*np.sin(theta_arr + zeta_spl(r)*np.sin(2*theta_arr)),
+ label='Miller')
+ax7.set_title('Flux surface shapes')
+ax7.set_xlabel('$R$/m')
+ax7.set_ylabel('$Z$/m')
+ax7.set_aspect('equal')
+ax7.legend(loc=10) #, prop={'size': 10})
+
+provide_conversions = 0
+if provide_conversions:
+ f = open('rho_tor-r_a_conv', 'w')
+ f.write('#r/a rho_tor\n')
+ for i in range(0, nw):
+ f.write('%16.8e %16.8e\n'%(ravg_spl(linpsi[i])/r_avg[-1], rho_tor_spl(linpsi[i])))
+
+if show_plots:
+ fig.tight_layout()
+ plt.show()
+
+# #Fourier series representation of flux surfaces
+# RZ_c=R+1j*Z
+# for i in range(0,nw,4):
+# RZ_c_fft=fft(RZ_c[i])
+# #plot(RZ_c_fft)
+# trunc_order=2
+# RZ_c_fft[trunc_order:nw/2]=0.;RZ_c_fft[nw/2+1:-trunc_order+1]=0.
+# RZ_trunc=ifft(RZ_c_fft)
+# R_trunc=real(RZ_trunc)
+# Z_trunc=imag(RZ_trunc)
+# plot(R_trunc,Z_trunc,lw=2,ls='--')
+# contour(rgrid,zgrid,psirz,20)
+# gca().set_aspect('equal')
+# show()
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index 09b285b87e7526eed4415cc0159c3430f000b7ea..215a920d5753933fa524a6cdb26377480ee81713 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -6,8 +6,8 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
default_plots_options
% Partition of the computer where the data have to be searched
% PARTITION='/Users/ahoffmann/gyacomo/results/paper_3/';
-% PARTITION='/misc/gyacomo23_outputs/paper_3/';
-PARTITION = '';
+PARTITION='/misc/gyacomo23_outputs/paper_3/';
+% PARTITION = '';
%% Paper 3
% resdir = 'DTT_rho85/3x2x192x48x32';
% resdir = 'DTT_rho85/3x2x192x48x32_NT';
@@ -15,8 +15,11 @@ PARTITION = '';
% resdir = 'DTT_rho98/3x2x192x48x32_0.25grad';
% resdir = 'LM_DIIID_rho95/5x3x512x92x32';
% resdir = 'LM_DIIID_rho95/3x2x512x92x32';
+% resdir = 'DIIID_LM_rho90/3x2x256x128x32';
+% resdir = 'DTT_rho85_geom_scan/P8_J4_delta_nuDGGK_conv_test/delta_-0.3_nu_0.9';
+resdir = 'NT_DIIID_Austin2019_rho95/3x2x256x64x32';
% resdir = '../testcases/cyclone_example';
-resdir = '../testcases/CBC_ExBshear/std';
+% resdir = '../testcases/CBC_ExBshear/std';
% resdir = '../results/paper_3/HM_DTT_rho98/3x2x128x64x64';
%%
J0 = 00; J1 = 10;
@@ -26,7 +29,7 @@ DATADIR = [PARTITION,resdir,'/'];
data = {};
data = compile_results_low_mem(data,DATADIR,J0,J1);
-if 0
+if 1
%% Plot transport and phi radial profile
[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
% [data.PSI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'psi');
@@ -60,7 +63,7 @@ options.NAME = 'N_i^{00}';
% options.NAME = 'Q_x';
% options.NAME = 'n_i';
% options.NAME = 'n_i-n_e';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
@@ -73,24 +76,25 @@ options.RESOLUTION = 256;
create_film(data,options,'.gif')
end
-if 1
+if 0
%% field snapshots
% Options
-[data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00');
-[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
+[data.Na00, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'Na00');
+[data.PHI, data.Ts3D] = compile_results_3D(data.folder,J0,J1,'phi');
data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
-options.INTERP = 1;
+options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
options.NORMALIZE = 0;
-options.LOGSCALE = 1;
+options.LOGSCALE = 0;
options.CLIMAUTO = 1;
options.NAME = 'N_i^{00}';
+% options.NAME = 's_{Ey}';
% options.NAME = '\phi';
-options.PLAN = 'kxky';
+options.PLAN = 'yz';
options.COMP = 'avg';
-options.TIME = [0 100];
+options.TIME = [50 100];
% options.TIME = data.Ts3D(1:2:end);
options.RESOLUTION = 256;
fig = photomaton(data,options);
@@ -127,12 +131,12 @@ options.TIME = data.Ts3D;
options.KX_TW = [ 0 20]; %kx Growth rate time window
options.KY_TW = [ 0 50]; %ky Growth rate time window
options.NMA = 1;
-options.NMODES = 3;
+options.NMODES = 50;
options.iz = 'avg'; % avg or index
options.ik = 9; % sum, max or index
options.fftz.flag = 0;
-options.FIELD = 'Ni00';
-% options.FIELD = 'phi';
+% options.FIELD = 'Ni00';
+options.FIELD = 'phi';
options.GOK2 = 0;
fig = mode_growth_meter(data,options);
% save_figure(data,fig,'.png')
diff --git a/wk/get_miller_GENE_py.m b/wk/get_miller_GENE_py.m
new file mode 100644
index 0000000000000000000000000000000000000000..056eae650cb9fb271e2282f3ff3ff9ad5deafd50
--- /dev/null
+++ b/wk/get_miller_GENE_py.m
@@ -0,0 +1,54 @@
+function [variables] = get_miller_GENE_py(prof_folder,rho)
+
+filePath = [prof_folder,'/equilibrium.txt'];
+
+command = ['python extract_miller_from_eqdsk.py ',...
+ filePath,' ',num2str(rho),' > tmp.txt'];
+system(command);
+
+% Specify the path to the output file
+filePath = 'tmp.txt';
+
+% Read the content of the file
+fileID = fopen(filePath, 'r');
+fileContent = textscan(fileID, '%s', 'Delimiter', '\n', 'Whitespace', '');
+fclose(fileID);
+fileContent = fileContent{1};
+
+% Initialize a structure to store variables
+variables = struct();
+
+% Loop through each line in the file
+for i = 1:length(fileContent)
+ line = strtrim(fileContent{i}); % Remove leading/trailing whitespaces
+
+ % Skip empty lines
+ if isempty(line)
+ continue;
+ end
+
+ % Split the line into variable name and value
+ parts = strsplit(line, '=');
+
+ % Skip lines that don't have '=' or have more than two parts
+ if length(parts) ~= 2
+ continue;
+ end
+
+ % Extract variable name and value
+ variableName = strtrim(parts{1});
+ variableValue = str2double(strtrim(parts{2}));
+
+ % Check if the conversion to double was successful
+ if isnan(variableValue)
+ % If conversion fails, try as a string
+ variableValue = strtrim(parts{2});
+ end
+
+ % Store the variable in the structure
+ variables.(genvarname(variableName)) = variableValue;
+end
+
+system('rm tmp.txt');
+
+end
\ No newline at end of file
diff --git a/wk/get_param_from_profiles.m b/wk/get_param_from_profiles.m
new file mode 100644
index 0000000000000000000000000000000000000000..6376d03511eb590356d6d43c633f33b129535e5b
--- /dev/null
+++ b/wk/get_param_from_profiles.m
@@ -0,0 +1,38 @@
+function [params, profiles] = get_param_from_profiles(folder,rho,Lref,mref,Bref,FROMPROFILE)
+%get_param_from_profiles compute the input param for GYACOMO from profiles
+%data
+m_e = 9.11e-31;
+params.SIGMA = sqrt(m_e/mref);
+if ~FROMPROFILE
+ nmvm = 5;
+ stencil = 9;
+ for field = {'ne','te','ti','wExB'}
+ f_.name = field{1};
+ fxy = load([folder,f_.name,'.txt']);
+ f_.x = movmean(fxy(:,1),nmvm);
+ f_.y = movmean(fxy(:,2),nmvm);
+ f_.K =-calculate_derivative(f_.x,f_.y,stencil)./f_.y;
+ profiles.(field{1}) = f_;
+ end
+ profiles.nmvm = nmvm;
+ profiles.stencil = stencil;
+else
+ profiles = read_DIIID_profile([folder,'/profiles.txt']);
+end
+
+ n_e = interpolate_at_x0(profiles.ne.x,profiles.ne.y,rho);
+ T_e = interpolate_at_x0(profiles.te.x,profiles.te.y,rho);
+ T_i = interpolate_at_x0(profiles.ti.x,profiles.ti.y,rho);
+ wExB = 0;%interpolate_at_x0(profiles.wExB.x,profiles.wExB.y,rho);
+ K_w = 0;%interpolate_at_x0(profiles.wExB.x,profiles.wExB.K,rho);
+ params.K_Ne = interpolate_at_x0(profiles.ne.x,profiles.ne.K,rho);
+ params.K_Ni = params.K_Ne;
+ params.K_Te = interpolate_at_x0(profiles.te.x,profiles.te.K,rho);
+ params.K_Ti = interpolate_at_x0(profiles.ti.x,profiles.ti.K,rho);
+ params.TAU = T_i/T_e;
+ params.nuGENE = 2.3031E-5*Lref*(n_e)/(T_e)^2*(24.-log(sqrt(n_e*1.0E13)/T_e*0.001));
+ params.NU = 3/8*sqrt(pi/2)*params.TAU.^(3/2)*params.nuGENE;
+ 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;
+end
\ No newline at end of file
diff --git a/wk/lin_run_script.m b/wk/lin_run_script.m
index 13da10b74bff69506a55c95ff094a89455ed646b..f872fc8767ddeeb77a2475406efafb4050956861 100644
--- a/wk/lin_run_script.m
+++ b/wk/lin_run_script.m
@@ -29,24 +29,51 @@ EXECNAME = 'gyacomo23_dp'; % double precision
% run lin_DIIID_LM_rho95
% run lin_JET_rho97
% run lin_Entropy
-run lin_ITG
+% run lin_ITG
% run lin_KBM
+% run lin_RHT
+rho = 0.95; TRIANG = 'NT'; READPROF = 0;
+prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/'];
+% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson/',TRIANG,'/'];
+% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/'];
+run lin_DIIID_data
+% run lin_STEP_EC_HD_psi71
+% run lin_STEP_EC_HD_psi49
+if 1
+% Plot the profiles
+ plot_params_vs_rho(geom,prof_folder,rho,0.5,1.1,Lref,mref,Bref,READPROF);
+end
%% Change parameters
% EXBRATE = 0.0; % Background ExB shear flow
-NU = 0.1;
-CO = 'SG';
+% K_Ti = 5.3;
+% NU = 0.001;
+CO = 'LD';
GKCO = 1;
-NY = 16;
-NX = 2;
-PMAX = 4;
+kymin= 0.3; LY = 2*pi/kymin;
+NY = 2;
+NX = 4;
+NZ = 32;
+PMAX = 2;
JMAX = PMAX/2;
-ky = 0.1; LY = 2*pi/ky;
-% DT = 1e-4;
-% TAU = 2.1;
-% SIGMA_E = 0.02;
-TMAX = 25;
+DT = 20e-4;
+EXBRATE = 0;
+% EPS = 0.25;
+% KAPPA = 1.0;
+S_DELTA = min(2.0,S_DELTA);
+% DELTA = -DELTA;
+% PT parameters
+% TAU = 5.45;
+% K_Ne = 0; %vs 2.79
+% K_Ni = 0;
+% K_Te = 9.6455;%vs 17.3
+% K_Ti = 3.3640;%vs 5.15
+% BETA = 7.9e-4;
+% NU = 0.1;
+MU_X = 0.0; MU_Y = 0.0;
+SIGMA_E = 0.04;
+TMAX = 1;
% DTSAVE0D = 200*DT;
-% DTSAVE3D = TMAX/50;
+DTSAVE3D = 0.1;
%%-------------------------------------------------------------------------
%% RUN
setup
@@ -55,8 +82,8 @@ setup
if RUN
MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
% RUN =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
- % RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
- RUN =['time ',mpirun,' -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
+ RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
+ % RUN =['time ',mpirun,' -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0;'];
% RUN =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;'];
% RUN =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
% RUN = ['./../../../bin/gyacomo23_sp 0;'];
@@ -75,16 +102,15 @@ data = {}; % Initialize data as an empty cell array
% load grids, inputs, and time traces
data = compile_results_low_mem(data,LOCALDIR,J0,J1);
-
-if 1 % Activate or not
+if 0 % Activate or not
%% plot mode evolution and growth rates
% Load phi
[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
options.NORMALIZED = 0;
options.TIME = data.Ts3D;
% Time window to measure the growth of kx/ky modes
-options.KX_TW = [0.2 1]*data.Ts3D(end);
-options.KY_TW = [0.2 1]*data.Ts3D(end);
+options.KY_TW = [0.5 1.0]*data.Ts3D(end);
+options.KX_TW = [0.5 1.0]*data.Ts3D(end);
options.NMA = 1; % Set NMA option to 1
options.NMODES = 999; % Set how much modes we study
options.iz = 'avg'; % Compressing z
@@ -92,20 +118,85 @@ options.ik = 1; %
options.GOK2 = 0; % plot gamma/k^2
options.fftz.flag = 0; % Set fftz.flag option to 0
options.FIELD = 'phi';
+options.SHOWFIG = 1;
[fig, kykx, wkykx, ekykx] = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
end
-if (1 && NZ>4)
+if (0 && NZ>4)
%% Ballooning plot
[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
if data.inputs.BETA > 0
[data.PSI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'psi');
end
options.time_2_plot = [120];
-options.kymodes = [0.2 0.3 0.4];
+options.kymodes = [100];
options.normalized = 1;
options.PLOT_KP = 0;
% options.field = 'phi';
-fig = plot_ballooning(data,options);
+options.SHOWFIG = 1;
+[fig, chi, phib, psib, ~] = plot_ballooning(data,options);
+end
+
+if 0
+%% RH TEST
+[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+ikx = 2; iky = 1; t0 = 1; t1 = data.Ts3D(end);
+[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));
+plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3));
+t_ = data.Ts3D(it0:it1);
+% TH = 1.635*EPS^1.5 + 0.5*EPS^2+0.35*EPS^2.5; theory = 1/(1+Q0^2*TH/EPS^2);
+clr_ = lines(20);
+figure
+plot(t_, plt(data.PHI)); hold on;
+plot(t_,0.5* exp(-t_*NU)+theory,'--k');
+plot([t_(1) t_(end)],theory*[1 1],'-k');
+plot([t_(1) t_(end)],mean(plt(data.PHI))*[1 1],'-g');
+xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
+title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.grids.kx(ikx),data.grids.ky(iky)))
+end
+
+if 0
+ %% Geometry
+ plot_metric(data);
end
+if 1
+ %% Compiled plot for lin EM analysis
+ [data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+ if data.inputs.BETA > 0
+ [data.PSI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'psi');
+ end
+ options.time_2_plot = [120];
+ options.kymodes = [100];
+ options.normalized = 1;
+ options.PLOT_KP = 0;
+ options.SHOWFIG = 0;
+ options.NORMALIZED = 0;
+ options.TIME = data.Ts3D;
+ % Time window to measure the growth of kx/ky modes
+ options.KX_TW = [0.5 1]*data.Ts3D(end);
+ options.KY_TW = [0.5 1]*data.Ts3D(end);
+ options.NMA = 1; % Set NMA option to 1
+ options.NMODES = 999; % Set how much modes we study
+ options.iz = 'avg'; % Compressing z
+ options.ik = 1; %
+ options.GOK2 = 0; % plot gamma/k^2
+ options.fftz.flag = 0; % Set fftz.flag option to 0
+ options.FIELD = 'phi';
+ options.SHOWFIG = 0;
+ [~, chi, phib, psib, ~] = plot_ballooning(data,options);
+ [~, kykx, wkykx, ekykx] = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
+ [~,~,R,Z] = plot_metric(data,options);
+ figure;
+ subplot(3,2,1); plot(chi,real(phib),'-b'); hold on;
+ plot(chi,imag(phib),'-r'); xlabel('$\chi$'); ylabel('$\phi$')
+ subplot(3,2,3); plot(chi,real(psib),'-b'); hold on;
+ plot(chi,imag(psib),'-r'); xlabel('$\chi$'); ylabel('$\psi$')
+ subplot(3,2,5); plot(squeeze(kykx(:,1)),squeeze(real(wkykx(:,1))),'o--'); hold on;
+ plot(squeeze(kykx(:,1)),squeeze(imag(wkykx(:,1))),'o--');
+ xlabel('$k_y\rho_s$'); ylabel('$\gamma,\omega$');
+ R = R*geom.R0; Z = Z*geom.R0;
+ subplot(1,2,2); plot(R,Z,'-k'); xlabel('$R$'); ylabel('$Z$');
+ axis equal
+ title(data.paramshort);
+end
\ No newline at end of file
diff --git a/wk/lin_scan_script.m b/wk/lin_scan_script.m
index 2dea07f623d52d7151946d1c7ed97270ef2d0c6c..bf36d44d2d156a2c46e53834f47edaa055c7b349 100644
--- a/wk/lin_scan_script.m
+++ b/wk/lin_scan_script.m
@@ -4,8 +4,8 @@
% for benchmarking and debugging purposes since it makes Matlab "busy".
%% Set up the paths for the necessary Matlab modules
-gyacomodir = pwd;
-gyacomodir = gyacomodir(1:end-2);
+wkdir = pwd;
+gyacomodir = wkdir(1:end-2);
mpirun = 'mpirun';
% mpirun = '/opt/homebrew/bin/mpirun'; % for macos
addpath(genpath([gyacomodir,'matlab'])) % Add matlab folder
@@ -18,26 +18,34 @@ addpath(genpath([gyacomodir,'wk/parameters'])) % Add parameters folder
RUN = 1; % To run or just to load
RERUN = 0; % rerun if the data does not exist
default_plots_options
-EXECNAME = 'gyacomo23_sp'; % single precision
-%EXECNAME = 'gyacomo23_dp'; % double precision
+% EXECNAME = 'gyacomo23_sp'; % single precision
+EXECNAME = 'gyacomo23_dp'; % double precision
%% Setup parameters
% run lin_DTT_AB_rho85
% run lin_DTT_AB_rho98
-run lin_JET_rho97
+% run lin_JET_rho97
% run lin_Entropy
% run lin_ITG
+% run lin_RHT
+rho = 0.95; TRIANG = 'NT'; READPROF = 1;
+% prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/'];
+% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson/',TRIANG,'/'];
+prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/'];
+run lin_DIIID_data
%% Change parameters
+NU = 1;
+TAU = 1;
NY = 2;
EXBRATE = 0;
-% SIGMA_E = 0.023;
+SIGMA_E = 0.023;
%% Scan parameters
SIMID = [SIMID,'_scan'];
-P_a = [2 4 6 8];
+P_a = [2 4];
% P_a = 2;
-ky_a = logspace(-1.5,1.5,30);
-CO = 'SG';
+ky_a = [0.01 0.02 0.05 0.1 0.2 0.5 1.0 2.0 5.0 10.0];
+CO = 'LD';
%% Scan loop
% arrays for the result
g_ky = zeros(numel(ky_a),numel(P_a));
@@ -50,10 +58,10 @@ for PMAX = P_a
i = 1;
for ky = ky_a
LY = 2*pi/ky;
- DT = 1e-5;%/(1+log(ky/0.05));%min(1e-2,1e-3/ky);
+ DT = 2e-4;%/(1+log(ky/0.05));%min(1e-2,1e-3/ky);
TMAX = 20;%min(10,1.5/ky);
DTSAVE0D = 0.1;
- DTSAVE3D = 0.1;
+ DTSAVE3D = 0.01;
%% RUN
setup
% naming
@@ -68,13 +76,13 @@ for PMAX = P_a
data_.outfilenames = [];
end
if RUN && (RERUN || isempty(data_.outfilenames) || (Ntime < 10))
- MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
+ MVIN =['cd ',LOCALDIR,';'];
% RUNG =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
RUNG =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
% RUNG =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;'];
% RUNG =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
% RUNG = ['./../../../bin/gyacomo23_sp 0;'];
- MVOUT='cd ../../../wk;';
+ MVOUT=['cd ',wkdir,';'];
system([MVIN,RUNG,MVOUT]);
end
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
@@ -95,11 +103,11 @@ 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 ...
[wkykx,ekykx] = compute_growth_rates(data_.PHI(:,:,:,it1:it2),data_.Ts3D(it1:it2));
- g_ky (i,j) = real(wkykx(2,1,end));
+ g_ky (i,j) = real(wkykx(2,1));
g_std(i,j) = real(ekykx(2,1));
- w_ky (i,j) = imag(wkykx(2,1,end));
+ w_ky (i,j) = imag(wkykx(2,1));
w_std(i,j) = imag(ekykx(2,1));
- [gmax, ikmax] = max(g_ky(i,j,:));
+ [gmax, ikmax] = max(g_ky(i,j));
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
end
diff --git a/wk/load_metadata_scan.m b/wk/load_metadata_scan.m
index 7ef57a678dfc60073e788d78e736c638f4d0a4be..a9df5c95be4a524979dbae0448bf4f1f2461c3c2 100644
--- a/wk/load_metadata_scan.m
+++ b/wk/load_metadata_scan.m
@@ -5,21 +5,13 @@ 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';
-
-% datafname = 'lin_KBM/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.03.mat';
-% datafname = 'lin_KBM/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.03.mat';
-% datafname = 'lin_AUG_scan/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.000152.mat';
-% datafname = 'lin_Entropy_scan/2x1_ky_0.05_0.75_P_2_8_DGDK_0_be_0.mat';
-% 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/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_8_kN_10_DGGK_0.1_be_0.0031.mat';
+% high density DIII-D
+% rho = 0.9
+% datafname = 'lin_DIIID_Oak_Nelson_high_density_PT_scan/6x32_ky_0.01_10_P_2_8_kN_0.32141_LDGK_0.0038027_be_0.0060039.mat';
+% datafname = 'lin_DIIID_Oak_Nelson_high_density_NT_scan/6x32_ky_0.01_10_P_2_8_kN_1.0883_LDGK_0.0080915_be_0.0015991.mat';
+% rho = 0.95
+% datagname = 'lin_DIIID_Oak_Nelson_high_density_PT_scan/6x32_ky_0.01_10_P_2_8_kN_0.62888_LDGK_0.0046858_be_0.0048708.mat';
+datafname = 'lin_DIIID_Oak_Nelson_high_density_NT_scan/6x32_ky_0.01_10_P_2_8_kN_1.6989_LDGK_0.0167_be_0.00075879.mat';
%% Chose if we filter gamma>0.05
FILTERGAMMA = 1;
diff --git a/wk/p3_geom_scan_analysis.m b/wk/p3_geom_scan_analysis.m
index 38d09ed95f1f5b1318aaa62ab72463775f9fdfbd..da097816abe751451d23eee47b426a038cda8c69 100644
--- a/wk/p3_geom_scan_analysis.m
+++ b/wk/p3_geom_scan_analysis.m
@@ -3,7 +3,7 @@ curdir = pwd;
partition= '/misc/gyacomo23_outputs/paper_3/';
% partition= '../results/paper_3/';
% Get the scan directory
-switch 1
+switch 2
case 1 % delta kappa scan
scandir = 'DTT_rho85_geom_scan/P2_J1_delta_kappa_scan/';
% scandir = 'DTT_rho85_geom_scan/P4_J2_delta_kappa_scan/';
@@ -15,21 +15,24 @@ switch 1
nml1 = 'GEOMETRY'; pnam1 = '$\hat s$'; attr1 = 'shear'; pref1 = 3.63;
nml2 = 'GEOMETRY'; pnam2 = '$q_0$'; attr2 = 'q0'; pref2 =-2.15;
case 3
- % scandir = 'DTT_rho85_geom_scan/P2_J1_delta_nuDGGK_scan/';
+ scandir = 'DTT_rho85_geom_scan/P2_J1_delta_nuDGGK_scan/';
% scandir = 'DTT_rho85_geom_scan/P4_J2_delta_nuDGGK_scan/';
+ % scandir = 'DTT_rho85_geom_scan/P8_J4_delta_nuDGGK_conv_test/';
% scandir = 'DTT_rho85_geom_scan/P2_J1_delta_nuSGGK_scan/';
% scandir = 'DTT_rho85_geom_scan/P4_J2_delta_nuSGGK_scan/';
+ % scandir = 'DTT_rho85_geom_scan/P8_J4_delta_nuSGGK_conv_test/';
+ % scandir = 'DTT_rho85_geom_scan/P4_J2_delta_nuSGGKii_scan/';
% scandir = 'DTT_rho85_geom_scan/P2_J1_delta_nuLDGK_scan/';
- scandir = 'DTT_rho85_geom_scan/P4_J2_delta_nuLDGK_scan/';
+ % scandir = 'DTT_rho85_geom_scan/P4_J2_delta_nuLDGK_scan/';
nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23;
- nml2 = 'MODEL'; pnam2 = '$\nu$'; attr2 = 'nu'; pref2 = 999;%0.5;
+ nml2 = 'MODEL'; pnam2 = '$\nu$'; attr2 = 'nu'; pref2 = 0.5;
end
scandir = [partition,scandir];
% Get a list of all items in the current directory
contents = dir(scandir);
% give ref value and param names
-REFVAL= 1;
+REFVAL= 0;
% normalize to the max
NORMAL= 0;
% Get and plot the fluxsurface
@@ -50,7 +53,7 @@ for i = 1:length(contents)
para2 = [para2 param.(nml2).(attr2)];
% Now you are in the subdirectory. You can perform operations here.
[t_all, Pxi_all, Qxi_all, Pxe_all, Qxe_all] = read_flux_out_XX(subdir);
- if(t_all(end) > 100)
+ if(t_all(end) > 50)
[fullAvg,sliceAvg,sliceErr] = sliceAverage(Qxi_all+Qxe_all,3);
Qxavg = [Qxavg fullAvg];
Qxerr = [Qxerr mean(sliceErr)];
@@ -126,7 +129,7 @@ if REFVAL
xref = XX(iref1,iref2);
yref = YY(iref1,iref2);
Qxref = Zavg(iref1,iref2);
- % Qrefname = ['$Q_{ref}=$',num2str(Qxref)];
+ Qrefname = ['$Q_{ref}=$',num2str(Qxref(1,1))];
else
Qxname = '$\langle Q_{tot} \rangle_t$';
end
@@ -148,9 +151,10 @@ else
imagesc_custom(xx_,yy_,Zavg'); hold on
CLIM = 'auto';
end
-% if REFVAL
-% plot(xref,yref,'xk','MarkerSize',14,'DisplayName',Qrefname)
-% end
+if REFVAL && ~((pref1==999) || (pref2==999))
+ plot(xref(1,1),yref(1,1),'xk','MarkerSize',14,'DisplayName',Qrefname)
+ legend('show')
+end
xlabel(pnam1); ylabel(pnam2);
title(Qxname)
colormap(bluewhitered); colorbar; clim(CLIM);
@@ -163,9 +167,9 @@ for i = 1:N2
'DisplayName',[pnam2,'=',num2str(para2(1,i))]);
hold on;
end
-% if REFVAL
-% plot(xref,0,'xk','MarkerSize',14,'DisplayName',Qrefname)
-% end
+if REFVAL && ~((pref1==999) || (pref2==999))
+ plot(xref(1,1),0,'xk','MarkerSize',14,'DisplayName',Qrefname)
+end
grid on
xlabel(pnam1); ylabel('$\langle Q_{tot} \rangle_t$');
legend('show','Location','northwest');
diff --git a/wk/parameters/lin_DIIID_AUSTIN.m b/wk/parameters/lin_DIIID_AUSTIN.m
new file mode 100644
index 0000000000000000000000000000000000000000..1ba1c97d776348099cb6ec8d802832fbe69d29f6
--- /dev/null
+++ b/wk/parameters/lin_DIIID_AUSTIN.m
@@ -0,0 +1,129 @@
+%% Reference values
+prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/'];
+% Get geometry parameters
+GEOMETRY= 'miller';
+if READPROF
+ geom = get_miller_GENE_py(prof_folder,rho);
+else
+ geom.q0 = 3.25; geom.shat = 3.65;
+ geom.kappa = 1.66; geom.s_kappa= 0.77;
+ geom.delta = 0.34; geom.s_delta= 1.13;
+ if strcmp(TRIANG,'NT')
+ geom.delta = -geom.delta;
+ end
+ geom.zeta =-0.038;geom.s_zeta =-0.13;
+ geom.trpeps = 1/3.15;
+ geom.Bref = 1;
+ geom.R0 = 3.15;
+ geom.a = 1;
+ geom.Lref = geom.R0;
+end
+Q0 = geom.q0; % safety factor
+SHEAR = geom.shat; % magnetic shear
+EPS = geom.trpeps; % inverse aspect ratio
+KAPPA = geom.kappa; % elongation
+S_KAPPA = geom.s_kappa;
+DELTA = geom.delta; % triangularity
+S_DELTA = geom.s_delta;
+ZETA = geom.zeta; % squareness
+S_ZETA = geom.s_zeta;
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+% See Austin et al. 2019, negative triangularity exp. in DIIID
+q_e= 1.602e-19;
+Z = 2;
+m_i= Z*1.67e-27;
+m_e= 9.11e-31;
+BT = geom.Bref; %T
+Ip = 0.9; %MA
+R0 = geom.R0; %m
+a0 = geom.a;%m
+P_in = 13; %input power in MW
+Lref = geom.Lref;
+Bref = abs(geom.Bref);
+mref = m_i;
+% Get values at a given location :
+[params,profiles] = get_param_from_profiles(prof_folder,rho,Lref,mref,Bref,READPROF);
+K_Ne = params.K_Ne;
+K_Ni = params.K_Ne;
+K_Te = params.K_Te;
+K_Ti = params.K_Ti;
+TAU = params.TAU;
+NU = params.NU;
+BETA = params.BETA;
+SIGMA_E= sqrt(m_e/mref);
+EXBRATE= params.EXBRATE;
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+PB_PHASE = 0;
+% Name of the simulation folder
+SIMID = ['lin_DIIID_Austin',TRIANG,'_',num2str(rho)];
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NA = 2; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+MHD_PD = 1;
+%% Set up grid parameters
+P = 2;
+J = P/2;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 6; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.3; % Size of the squared frequency domain in y direction
+NZ = 32; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+%% TIME PARAMETERS
+TMAX = 15; % Maximal time unit
+DT = 1e-4; % Time step
+DTSAVE0D = 0.5; % Sampling time for 0D arrays
+DTSAVE2D = -1; % Sampling time for 2D arrays
+DTSAVE3D = 0.5; % Sampling time for 3D arrays
+DTSAVE5D = TMAX/2; % Sampling time for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 1; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 5.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 1.0e-5; % Initial noise amplitude
+BCKGD0 = 0.0e-5; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1; % Cutoff for collision operator
+ADIAB_I = 0; % adiabatic ion model
\ No newline at end of file
diff --git a/wk/parameters/lin_DIIID_data.m b/wk/parameters/lin_DIIID_data.m
new file mode 100644
index 0000000000000000000000000000000000000000..297f11e3c8565ea031fff4d18e0bcd3850fc0f0e
--- /dev/null
+++ b/wk/parameters/lin_DIIID_data.m
@@ -0,0 +1,132 @@
+%% Reference values
+% prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/'];
+% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_1/',TRIANG,'/'];
+% Get geometry parameters
+GEOMETRY= 'miller';
+if READPROF
+ geom = get_miller_GENE_py(prof_folder,rho);
+else
+ % Ad hoc geometry parameters
+ geom.q0 = 4.8; geom.shat = 2.55;
+ geom.kappa = 1.57; geom.s_kappa= 0.48;
+ geom.delta = 0.4; geom.s_delta= 0.25;
+ if strcmp(TRIANG,'NT')
+ geom.delta = -geom.delta;
+ geom.s_delta = -geom.s_delta;
+ end
+ geom.zeta = 0.0; geom.s_zeta = 0.0;
+ geom.Bref = 1;
+ geom.R0 = 3.15;
+ geom.trpeps = rho/geom.R0;
+ geom.a = 1;
+ geom.Lref = geom.R0;
+end
+Q0 = geom.q0; % safety factor
+SHEAR = geom.shat; % magnetic shear
+EPS = geom.trpeps; % inverse aspect ratio
+KAPPA = geom.kappa; % elongation
+S_KAPPA = geom.s_kappa;
+DELTA = geom.delta; % triangularity
+S_DELTA = min(2.0,geom.s_delta);
+ZETA = geom.zeta; % squareness
+S_ZETA = geom.s_zeta;
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+% See Austin et al. 2019, negative triangularity exp. in DIIID
+q_e= 1.602e-19;
+Z = 1;
+m_i= 2*1.67e-27; % Deuterium
+m_e= 9.11e-31;
+BT = geom.Bref; %T
+Ip = 0.9; %MA
+R0 = geom.R0; %m
+a0 = geom.a;%m
+P_in = 13; %input power in MW
+Lref = geom.Lref;
+Bref = abs(geom.Bref);
+mref = m_i;
+% Get values at a given location :
+[params,profiles] = get_param_from_profiles(prof_folder,rho,Lref,mref,Bref,READPROF);
+K_Ne = params.K_Ne;
+K_Ni = params.K_Ne;
+K_Te = params.K_Te;
+K_Ti = params.K_Ti;
+TAU = params.TAU;
+BETA = params.BETA;
+SIGMA_E= sqrt(m_e/mref);
+NU = 3/8*sqrt(pi/2)*TAU^(3/2)*params.nuGENE;
+EXBRATE= params.EXBRATE;
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+PB_PHASE = 0;
+% Name of the simulation folder
+SIMID = ['lin_',prof_folder(21:end-4),'_',TRIANG];
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NA = 2; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+MHD_PD = 1;
+%% Set up grid parameters
+P = 2;
+J = P/2;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 6; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.3; % Size of the squared frequency domain in y direction
+NZ = 32; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+%% TIME PARAMETERS
+TMAX = 15; % Maximal time unit
+DT = 1e-4; % Time step
+DTSAVE0D = 0.5; % Sampling time for 0D arrays
+DTSAVE2D = -1; % Sampling time for 2D arrays
+DTSAVE3D = 0.5; % Sampling time for 3D arrays
+DTSAVE5D = TMAX/2; % Sampling time for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 1; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 5.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 1.0e-5; % Initial noise amplitude
+BCKGD0 = 0.0e-5; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1; % Cutoff for collision operator
+ADIAB_I = 0; % adiabatic ion model
\ No newline at end of file
diff --git a/wk/parameters/lin_STEP_EC_HD_psi49.m b/wk/parameters/lin_STEP_EC_HD_psi49.m
new file mode 100644
index 0000000000000000000000000000000000000000..d1fa57efe1e0dd57271a8282f9f1d03ea412b7f0
--- /dev/null
+++ b/wk/parameters/lin_STEP_EC_HD_psi49.m
@@ -0,0 +1,116 @@
+% Taken from Table 2, Kennedy et al. 2023 preprint on STEP linear
+% simulations. rho = 0.64 (psi=0.49)
+%% Reference values
+% Get geometry parameters
+GEOMETRY= 'miller';
+Q0 = 3.6; % safety factor
+SHEAR = 1.20; % magnetic shear
+EPS = 0.56; % inverse aspect ratio
+KAPPA = 2.56; % elongation
+S_KAPPA = 0.06;
+DELTA = 0.29; % triangularity
+S_DELTA = 0.46;
+ZETA = 0; % squareness
+S_ZETA = 0;
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+% See Austin et al. 2019, negative triangularity exp. in DIIID
+q_e= 1.602e-19;
+Z = 2;
+m_i= Z*1.67e-27;
+m_e= 9.11e-31;
+BT = 3.2; %T
+Ip = 20.9; %MA
+R0 = 3.6; %m
+a0 = 2;%m
+P_in = 300; %input power in MW
+Lref = R0;
+Bref = abs(BT);
+mref = m_i;
+T_e = 18.0; % KeV
+n_e = 20.05; %10^19 m-3
+% Get values at a given location :
+K_Ne = 1.06;
+K_Ni = K_Ne;
+K_Te = 1.58;
+K_Ti = 1.82;
+TAU = 1.0;
+nuGENE = 2.3031E-5*Lref*(n_e)/(T_e)^2*(24.-log(sqrt(n_e*1.0E13)/T_e*0.001));
+NU = 0.05;%3/8*sqrt(pi/2)*sqrt(TAU)*nuGENE;
+BETA = 0.09;%403.0E-5*n_e*T_e/(Bref*Bref);
+SIGMA_E= sqrt(m_e/mref);
+EXBRATE= 0;
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+PB_PHASE = 0;
+% Name of the simulation folder
+SIMID = 'lin_STEP_EC_HD_psi49';
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NA = 2; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+MHD_PD = 1;
+%% Set up grid parameters
+P = 2;
+J = P/2;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 6; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.3; % Size of the squared frequency domain in y direction
+NZ = 32; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+%% TIME PARAMETERS
+TMAX = 15; % Maximal time unit
+DT = 1e-4; % Time step
+DTSAVE0D = 0.5; % Sampling time for 0D arrays
+DTSAVE2D = -1; % Sampling time for 2D arrays
+DTSAVE3D = 0.5; % Sampling time for 3D arrays
+DTSAVE5D = 100; % Sampling time for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 1; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 5.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 1.0e-5; % Initial noise amplitude
+BCKGD0 = 0.0e-5; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1; % Cutoff for collision operator
+ADIAB_I = 0; % adiabatic ion model
\ No newline at end of file
diff --git a/wk/parameters/lin_STEP_EC_HD_psi71.m b/wk/parameters/lin_STEP_EC_HD_psi71.m
new file mode 100644
index 0000000000000000000000000000000000000000..9881468e8407551ba05f320b92e361239507b99a
--- /dev/null
+++ b/wk/parameters/lin_STEP_EC_HD_psi71.m
@@ -0,0 +1,116 @@
+% Taken from Table 2, Kennedy et al. 2023 preprint on STEP linear
+% simulations. rho = 0.7
+%% Reference values
+% Get geometry parameters
+GEOMETRY= 'miller';
+Q0 = 4.0; % safety factor
+SHEAR = 1.56; % magnetic shear
+EPS = 0.56; % inverse aspect ratio
+KAPPA = 2.57; % elongation
+S_KAPPA = 0.19;
+DELTA = 0.32; % triangularity
+S_DELTA = 0.54;
+ZETA = 0; % squareness
+S_ZETA = 0;
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+% See Austin et al. 2019, negative triangularity exp. in DIIID
+q_e= 1.602e-19;
+Z = 2;
+m_i= Z*1.67e-27;
+m_e= 9.11e-31;
+BT = 3.2; %T
+Ip = 20.9; %MA
+R0 = 3.6; %m
+a0 = 2;%m
+P_in = 300; %input power in MW
+Lref = R0;
+Bref = abs(BT);
+mref = m_i;
+T_e = 18.0; % KeV
+n_e = 20.05; %10^19 m-3
+% Get values at a given location :
+K_Ne = 1.54;
+K_Ni = K_Ne;
+K_Te = 1.77;
+K_Ti = 1.96;
+TAU = 1.0;
+nuGENE = 2.3031E-5*Lref*(n_e)/(T_e)^2*(24.-log(sqrt(n_e*1.0E13)/T_e*0.001));
+NU = 0.05;%3/8*sqrt(pi/2)*sqrt(TAU)*nuGENE;
+BETA = 0.07;%403.0E-5*n_e*T_e/(Bref*Bref);
+SIGMA_E= sqrt(m_e/mref);
+EXBRATE= 0;
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+PB_PHASE = 0;
+% Name of the simulation folder
+SIMID = 'lin_STEP_EC_HD_psi71';
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NA = 2; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+MHD_PD = 1;
+%% Set up grid parameters
+P = 2;
+J = P/2;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 6; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.3; % Size of the squared frequency domain in y direction
+NZ = 32; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+%% TIME PARAMETERS
+TMAX = 15; % Maximal time unit
+DT = 1e-4; % Time step
+DTSAVE0D = 0.5; % Sampling time for 0D arrays
+DTSAVE2D = -1; % Sampling time for 2D arrays
+DTSAVE3D = 0.5; % Sampling time for 3D arrays
+DTSAVE5D = 100; % Sampling time for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 1; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 5.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 1.0e-5; % Initial noise amplitude
+BCKGD0 = 0.0e-5; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1; % Cutoff for collision operator
+ADIAB_I = 0; % adiabatic ion model
\ No newline at end of file
diff --git a/wk/plot_params_vs_rho.m b/wk/plot_params_vs_rho.m
new file mode 100644
index 0000000000000000000000000000000000000000..0b3d64860fbba428210049054fa658e9b577785f
--- /dev/null
+++ b/wk/plot_params_vs_rho.m
@@ -0,0 +1,97 @@
+function [ ] = plot_params_vs_rho(geom,prof_folder,rho,rho_min,rho_max,Lref,mref,Bref,FROMPROFILE)
+
+if FROMPROFILE
+ % profiles = read_DIIID_profile([prof_folder,'/profiles.txt']);
+ [params,profiles] = get_param_from_profiles(prof_folder,rho,Lref,mref,Bref,FROMPROFILE);
+ m_e = 9.11e-31;
+ sigma = sqrt(m_e/mref);
+ TAU_a = profiles.ti.y./profiles.te.y;
+ nuGENE = 2.3031E-5*Lref*(profiles.ne.y)./(profiles.te.y).^2 ...
+ .*(24.-log(sqrt(profiles.ne.y*1.0E13)./profiles.te.y*0.001));
+ NU_a = 3/8*sqrt(pi/2)*TAU_a.^(3/2).*nuGENE;
+ BETA_a = 403.0E-5*profiles.ne.y.*profiles.te.y/(Bref*Bref);
+ figure
+ subplot(231)
+ plot(profiles.ne.x,profiles.ne.y); hold on
+ xlabel('$\rho$'); ylabel('$n_e(\times 10^{19}m^{-1})$')
+ xlim([rho_min rho_max]); xline(rho,'--k','Linewidth',2,'DisplayName','$\rho_{sim}$')
+ grid on
+ subplot(234)
+ plot(profiles.te.x,profiles.te.y,'DisplayName','electrons'); hold on
+ plot(profiles.ti.x,profiles.ti.y,'DisplayName','ions'); hold on
+ xlabel('$\rho$'); ylabel('$T$(keV)')
+ xlim([rho_min rho_max]); xline(rho,'--k','Linewidth',2,'DisplayName','$\rho_{sim}$')
+ grid on
+ subplot(235)
+ yyaxis("left")
+ semilogy(profiles.ti.x,NU_a); hold on
+ ylabel('$\nu$')
+ yyaxis("right")
+ plot(profiles.ti.x,BETA_a*100); hold on
+ ylabel('$\beta$[\%]')
+ xlim([rho_min rho_max]);xlabel('$\rho$'); xline(rho,'--k','Linewidth',2,'DisplayName','$\rho_{sim}$')
+ grid on
+ subplot(232)
+ plot(profiles.te.x,profiles.te.K,'DisplayName','$\kappa_{Te}$'); hold on
+ plot(profiles.ti.x,profiles.ti.K,'DisplayName','$\kappa_{Ti}$'); hold on
+ plot(profiles.ne.x,profiles.ne.K,'DisplayName','$\kappa_{n}$'); hold on
+ xlim([rho_min rho_max]); xline(rho,'--k','Linewidth',2,'DisplayName','$\rho_{sim}$')
+ xlabel('$\rho$'); ylabel('$\kappa_\chi=R_0/L_\chi$')
+ legend show
+ grid on
+else
+ rho_a = linspace(rho_min,rho_max,100);
+
+ NU_a = zeros(size(rho_a));
+ BETA_a = zeros(size(rho_a));
+ K_Ne_a = zeros(size(rho_a));
+ K_Ti_a = zeros(size(rho_a));
+ K_Te_a = zeros(size(rho_a));
+
+ for i = 1:numel(rho_a)
+ rho_ = rho_a(i);
+ [params,profiles] = get_param_from_profiles(prof_folder,rho_,Lref,mref,Bref,FROMPROFILE);
+ NU_a(i) = 3/8*sqrt(pi/2)*params.TAU.^(3/2)*params.nuGENE;
+ BETA_a(i) = params.BETA;
+ K_Ne_a(i) = params.K_Ne;
+ K_Ti_a(i) = params.K_Ti;
+ K_Te_a(i) = params.K_Te;
+ end
+ figure
+ subplot(231)
+ plot(profiles.ne.x,profiles.ne.y); hold on
+ xlabel('$\rho$'); ylabel('$n_e(\times 10^{19}m^{-1})$')
+ xlim([rho_min rho_max]); xline(rho,'--k','Linewidth',2,'DisplayName','$\rho_{sim}$')
+ grid on
+ subplot(234)
+ plot(profiles.te.x,profiles.te.y,'DisplayName','electrons'); hold on
+ plot(profiles.ti.x,profiles.ti.y,'DisplayName','ions'); hold on
+ xlabel('$\rho$'); ylabel('$T$(keV)')
+ xlim([rho_min rho_max]); xline(rho,'--k','Linewidth',2,'DisplayName','$\rho_{sim}$')
+ grid on
+ subplot(235)
+ yyaxis("left")
+ semilogy(rho_a,NU_a); hold on
+ ylabel('$\nu$')
+ yyaxis("right")
+ plot(rho_a,BETA_a*100); hold on
+ ylabel('$\beta$[\%]')
+ xlim([rho_min rho_max]);xlabel('$\rho$'); xline(rho,'--k','Linewidth',2,'DisplayName','$\rho_{sim}$')
+ grid on
+ subplot(232)
+ plot(rho_a,K_Te_a,'DisplayName','$\kappa_{Te}$'); hold on
+ plot(rho_a,K_Ti_a,'DisplayName','$\kappa_{Ti}$'); hold on
+ plot(rho_a,K_Ne_a,'DisplayName','$\kappa_{n}$'); hold on
+ xlim([rho_min rho_max]); xline(rho,'--k','Linewidth',2,'DisplayName','$\rho_{sim}$')
+ xlabel('$\rho$'); ylabel('$\kappa_\chi=R_0/L_\chi$')
+ legend show
+ grid on
+end
+ subplot(133)
+ [R,Z] = plot_miller(geom,rho,32,0);
+ plot(R,Z,'-b');
+ xlabel('R [m]');
+ ylabel('Z [m]');
+ axis tight
+ axis equal
+end
\ No newline at end of file
diff --git a/wk/read_DIIID_profile.m b/wk/read_DIIID_profile.m
new file mode 100644
index 0000000000000000000000000000000000000000..fbe1d617a89c8ae792d94fa70d3283b04c428984
--- /dev/null
+++ b/wk/read_DIIID_profile.m
@@ -0,0 +1,48 @@
+function data = read_DIIID_profile(filePath)
+ if exist([filePath,'.mat'])
+ tmp = load([filePath,'.mat']);
+ data = tmp.data;
+ else
+ % Read the content of the file
+ fileID = fopen(filePath, 'r');
+ fileContent = textscan(fileID, '%s', 'Delimiter', '\n', 'Whitespace', '');
+ fclose(fileID);
+ fileContent = fileContent{1};
+
+ % Initialize variables
+ data = struct();
+
+ % Loop through each line in the file
+ for i = 1:length(fileContent)
+ line = strtrim(fileContent{i}); % Remove leading/trailing whitespaces
+
+ % Check if the line starts with '201'
+ if startsWith(line, '201')
+ % Extract column titles with potential units
+ titles = strsplit(line, ' ');
+ titles(cellfun(@isempty, titles)) = []; % Remove empty cells
+ columnTitle = titles{3}; % Assuming the data starts from the 3rd column
+
+ % Remove potential units from the column title
+ cleanTitle = regexprep(columnTitle, '\(\S*\)', ''); % Remove everything inside parentheses
+
+ % Initialize data field for the current variable
+ data.(genvarname(cleanTitle)).x = [];
+ data.(genvarname(cleanTitle)).y = [];
+ data.(genvarname(cleanTitle)).K = [];
+
+ % Extract data for the current column
+ data.(genvarname(cleanTitle)).psinorm = sscanf(line, '%f %f %f', [3, inf])';
+ else
+ % Extract data for the current section
+ currentSectionData = sscanf(line, '%f %f %f');
+ data.(genvarname(cleanTitle)).x = [data.(genvarname(cleanTitle)).x; currentSectionData(1)];
+ data.(genvarname(cleanTitle)).y = [data.(genvarname(cleanTitle)).y; currentSectionData(2)];
+ data.(genvarname(cleanTitle)).K = [data.(genvarname(cleanTitle)).K;-currentSectionData(3)];
+ end
+ end
+ % unit conversions
+ data.ne.y = data.ne.y * 10; % convert to 10^19 m-3
+ save([filePath,'.mat'], 'data');
+ end
+end