diff --git a/matlab/compute/compute_fluxtube_growth_rate.m b/matlab/compute/compute_fluxtube_growth_rate.m
index 4cae40dab4361d3f0fd7679a794dbbbd4cd8aef1..d0beda9da004f60fa15035417bb9a6db523e1c74 100644
--- a/matlab/compute/compute_fluxtube_growth_rate.m
+++ b/matlab/compute/compute_fluxtube_growth_rate.m
@@ -80,9 +80,16 @@ end
if PLOT > 2
xlabel([]); xticks([]);
subplot(2,2,4)
- semilogy(DATA.Ts3D,squeeze(abs(DATA.PHI(2,1,DATA.Nz/2,:)))); hold on;
- xlabel('t'); ylabel('$|\phi_{ky}|(t)$')
-
+ [~,ikx0] = min(abs(DATA.kx));
+ for i_ = 1:(DATA.Nkx+1)/2
+ iky = 1 + (DATA.ky(1) == 0);
+ ikx = ikx0 + (i_-1);
+ semilogy(DATA.Ts3D,squeeze(abs(DATA.PHI(iky,ikx,DATA.Nz/2,:))),...
+ 'DisplayName',['$k_x=',num2str(DATA.kx(ikx)),'$, $k_y=',num2str(DATA.ky(iky)),'$']);
+ hold on;
+ xlabel('t,'); ylabel('$|\phi_{ky}|(t)$')
+ end
+legend('show')
end
end
\ No newline at end of file
diff --git a/matlab/compute/ifourier_GENE.m b/matlab/compute/ifourier_GENE.m
index 3c3a454a2032a384d73508347e8301d7f618454e..1538c36dfb1cde9ae4803a82e3ebf32b7ecddaea 100644
--- a/matlab/compute/ifourier_GENE.m
+++ b/matlab/compute/ifourier_GENE.m
@@ -11,19 +11,19 @@ ny=2*nky-1;
if ny~=1
%note, we need one extra point which we set to zero for the ifft
- spectrumKxKyZ=zeros(ny,nx,nz);
- spectrumKxKyZ(1:nky,:,:)=field_c(:,:,:);
- spectrumKxKyZ((nky+1):(ny),1,:)=conj(field_c(nky:-1:2,1,:));
- spectrumKxKyZ((nky+1):(ny),2:nx,:)=conj(field_c(nky:-1:2,nx:-1:2,:));
+ spectrumKyKxZ=zeros(ny,nx,nz);
+ spectrumKyKxZ(1:nky,:,:)=field_c(:,:,:);
+ spectrumKyKxZ((nky+1):(ny),1,:)=conj(field_c(nky:-1:2,1,:));
+ spectrumKyKxZ((nky+1):(ny),2:nx,:)=conj(field_c(nky:-1:2,nx:-1:2,:));
else
%pad with zeros to interpolate on fine scale
ny=20;
- spectrumKxKyZ=zeros(nx,ny,nz);
- spectrumKxKyZ(:,2,:)=field_c(:,:,:);
+ spectrumKyKxZ=zeros(nx,ny,nz);
+ spectrumKyKxZ(:,2,:)=field_c(:,:,:);
end
%inverse fft, symmetric as we are using real data
-spectrumXKyZ=nx*ifft(spectrumKxKyZ,[],1);
+spectrumXKyZ=nx*ifft(spectrumKyKxZ,[],1);
field_r=ny*ifft(spectrumXKyZ,[],2,'symmetric');
clear spectrumKxKyZ
diff --git a/matlab/compute/process.m b/matlab/compute/process.m
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/matlab/dbg_zBC_map.m b/matlab/dbg_zBC_map.m
new file mode 100644
index 0000000000000000000000000000000000000000..48b1bffa6ce835e64293fa6ac9ea091502c1bced
--- /dev/null
+++ b/matlab/dbg_zBC_map.m
@@ -0,0 +1,70 @@
+Nkx = 16;
+Nky = 2;
+my = 0:(Nky-1);
+mx = zeros(1,Nkx);
+
+for ix = 1:Nkx
+ if(mod(Nkx,2) == 0)%even
+ mx_max = (Nkx/2);
+ if(ix<=Nkx/2+1)
+ mx(ix) = (ix-1);
+ else
+ mx(ix) = ix-Nkx-1;
+ end
+ else %odd
+ mx_max = (Nkx-1)/2;
+ if(ix<=(Nkx-1)/2+1)
+ mx(ix) = (ix-1);
+ else
+ mx(ix) = ix-Nkx-1;
+ end
+ end
+end
+disp(mx)
+
+Npol = 1;
+Nexc = 1;
+shear = 0.8;
+Ly = 120;
+dky = 2*pi/Ly;
+dkx = 2*pi*shear*dky/Nexc;
+
+kx = mx*dkx;
+ky = my*dky;
+
+kx_max = mx_max*dkx;
+ikx_zBC_R = zeros(Nky,Nkx);
+for iy = 1:Nky
+ shift = 2*pi*shear*ky(iy)*Npol;
+ for ix = 1:Nkx
+ kx_shift = kx(ix) + shift;
+ if 0%(kx_shift > kx_max)
+ ikx_zBC_R(iy,ix) = nan;
+ else
+ ikx_zBC_R(iy,ix) = ix+(iy-1)*Nexc;
+ if(ikx_zBC_R(iy,ix) > Nkx)
+ ikx_zBC_R(iy,ix) = ikx_zBC_R(iy,ix) - Nkx;
+ end
+ end
+ end
+end
+disp(ikx_zBC_R)
+
+kx_min = (-mx_max+(1-mod(Nkx,2)))*dkx;
+ikx_zBC_L = zeros(Nky,Nkx);
+for iy = 1:Nky
+ shift = 2*pi*shear*ky(iy)*Npol;
+ for ix = 1:Nkx
+ kx_shift = kx(ix) - shift;
+ if(kx_shift < kx_min)
+ ikx_zBC_L(iy,ix) = nan;
+ else
+ ikx_zBC_L(iy,ix) = ix-(iy-1)*Nexc;
+ if(ikx_zBC_L(iy,ix) < 1)
+ ikx_zBC_L(iy,ix) = ikx_zBC_L(iy,ix) + Nkx;
+ end
+ end
+ end
+end
+disp(ikx_zBC_L)
+
diff --git a/matlab/load/load_params.m b/matlab/load/load_params.m
index 720e243fc9e6921cdf12a30a14c1fa41952d5855..d3d141fbbe3e4a14d028505805cd34420b86347a 100644
--- a/matlab/load/load_params.m
+++ b/matlab/load/load_params.m
@@ -9,12 +9,20 @@ end
try
DATA.K_N = h5readatt(filename,'/data/input','K_n');
catch
- DATA.K_N = h5readatt(filename,'/data/input','k_N');
+ try
+ DATA.K_N = h5readatt(filename,'/data/input','k_N');
+ catch
+ DATA.K_N = h5readatt(filename,'/data/input','k_Ni');
+ end
end
try
DATA.K_T = h5readatt(filename,'/data/input','K_T');
catch
- DATA.K_T = h5readatt(filename,'/data/input','k_T');
+ try
+ DATA.K_T = h5readatt(filename,'/data/input','k_T');
+ catch
+ DATA.K_T = h5readatt(filename,'/data/input','k_Ti');
+ end
end
DATA.Q0 = h5readatt(filename,'/data/input','q0');
DATA.SHEAR = h5readatt(filename,'/data/input','shear');
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index 6e6df0128c30bf11d8da6142d221ca877a42a098..a37b3b714f623306f9193552b446ca552916f681 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -13,8 +13,12 @@ function [FIG] = plot_ballooning(data,options)
psi_imag=imag(data.PSI(:,:,:,it1));
ncol = 2;
end
- % Apply baollooning tranform
- nexc = round(data.ky(2)*data.SHEAR*2*pi/data.kx(2));
+ % Apply ballooning transform
+ if(data.Nkx > 1)
+ nexc = round(data.ky(2)*data.SHEAR*2*pi/data.kx(2));
+ else
+ nexc = 1;
+ end
for iky=ikyarray
dims = size(phi_real);
Nkx = dims(2);
@@ -22,12 +26,18 @@ function [FIG] = plot_ballooning(data,options)
Npi = (Nkx-1)-2*nexc*(is-1);
if(Npi <= 1)
ordered_ikx = 1;
- else
+ elseif(mod(Nkx,2) == 0)
tmp_ = (Nkx-is+1):-is:(Nkx/2+2);
ordered_ikx = [tmp_(end:-1:1), 1:is:(Nkx/2)];
+ else
+ Np_ = (Nkx+1)/(2*is);
+ ordered_ikx = [(Np_+1):Nkx 1:Np_];
+ end
+ try
+ idx=data.kx./data.kx(2);
+ catch
+ idx=0;
end
-
- idx=data.kx./data.kx(2);
idx= idx(ordered_ikx);
Nkx = numel(idx);
@@ -36,7 +46,7 @@ function [FIG] = plot_ballooning(data,options)
b_angle = phib_real;
for i_ =1:Nkx
- start_ = (i_ -1)*dims(3) +1;
+ start_ = (i_-1)*dims(3) +1;
end_ = i_*dims(3);
ikx = ordered_ikx(i_);
phib_real(start_:end_) = phi_real(iky,ikx,:);
@@ -47,7 +57,7 @@ function [FIG] = plot_ballooning(data,options)
coordz = data.z;
for i_ =1: Nkx
for iz=1:dims(3)
- ii = dims(3)*(i_ -1) + iz;
+ ii = dims(3)*(i_-1) + iz;
b_angle(ii) =coordz(iz) + 2*pi*idx(i_)./is;
end
end
@@ -56,7 +66,7 @@ function [FIG] = plot_ballooning(data,options)
% normalize real and imaginary parts at chi =0
if options.normalized
[~,idxLFS] = min(abs(b_angle -0));
- normalization = abs(phib( idxLFS));
+ normalization = (phib( idxLFS));
else
normalization = 1;
end
@@ -73,8 +83,7 @@ function [FIG] = plot_ballooning(data,options)
title(['$k_y=',sprintf('%2.2f',data.ky(iky)),...
',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']);
- if data.BETA > 0
-
+ if data.BETA > 0
psib_real = zeros( Nkx*dims(3) ,1);
psib_imag = psib_real;
for i_ =1:Nkx
@@ -84,11 +93,14 @@ function [FIG] = plot_ballooning(data,options)
psib_real(start_:end_) = psi_real(iky,ikx,:);
psib_imag(start_:end_) = psi_imag(iky,ikx,:);
end
-
+ psib = psib_real(:) + 1i * psib_imag(:);
+ psib_norm = psib / normalization;
+ psib_real_norm = real( psib_norm);
+ psib_imag_norm = imag( psib_norm);
subplot(numel(ikyarray),ncol,ncol*(iplot-1)+2)
- plot(b_angle/pi, psib_real/ normalization,'-b'); hold on;
- plot(b_angle/pi, psib_imag/ normalization ,'-r');
- plot(b_angle/pi, sqrt(psib_real .^2 + psib_imag.^2)/ normalization,'-k');
+ 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')
xlabel('$\chi / \pi$')
ylabel('$\psi_B (\chi)$');
diff --git a/matlab/plot/plot_gbms_ballooning.m b/matlab/plot/plot_gbms_ballooning.m
new file mode 100644
index 0000000000000000000000000000000000000000..c995f3c8128b9e071d8afe896b0e47703213e1e2
--- /dev/null
+++ b/matlab/plot/plot_gbms_ballooning.m
@@ -0,0 +1,74 @@
+function [ ] = plot_gbms_ballooning(resfile)
+% perform ballooning transformation from phi.dat.h5 file
+
+% read data and attributes
+coordkx = h5read(resfile,'/data/var2d/phi/coordkx');
+Nkx = (length(coordkx)-1)/2;
+coordz = h5read(resfile,'/data/var2d/phi/coordz');
+Nz = length(coordz);
+coordtime =h5read(resfile,'/data/var2d/time');
+Nt = length(coordtime) ;
+
+iframe = Nt -1;
+dataset = ['/data/var2d/phi/',num2str(iframe,'%06d')];
+phi = h5read(resfile,dataset);
+try
+dataset = ['/data/var2d/psi/',num2str(iframe,'%06d')];
+psi = h5read(resfile,dataset);
+catch
+ psi = 0;
+end
+% Apply baollooning tranform
+dims = size(phi.real);
+phib_real = zeros( dims(1)*Nz ,1);
+phib_imag= phib_real;
+psib_real= phib_real;
+psib_imag= phib_real;
+b_angle = phib_real;
+
+midpoint = floor((dims(1)*Nz )/2)+1;
+
+for ip =1: dims(1)
+ start_ = (ip -1)*Nz +1;
+ end_ = ip*Nz;
+ phib_real(start_:end_) = phi.real(ip,:);
+ phib_imag(start_:end_) = phi.imaginary(ip,:);
+ try
+ psib_real(start_:end_) = psi.real(ip,:);
+ psib_imag(start_:end_) = psi.imaginary(ip,:);
+ catch
+ psib_real(start_:end_) = 0;
+ psib_imag(start_:end_) = 0;
+ end
+end
+
+% Define ballooning angle
+idx = -Nkx:1:Nkx;
+for ip =1: dims(1)
+ for iz=1:Nz
+ ii = Nz*(ip -1) + iz;
+ b_angle(ii) =coordz(iz) + 2*pi*idx(ip);
+ end
+end
+
+% normalize real and imaginary parts at chi =0
+[~,idxLFS] = min(abs(b_angle -0));
+phib = phib_real + 1i*phib_imag;
+psib = psib_real + 1i*psib_imag;
+% normalize to the outer mid-plane
+norm = (phib(idxLFS));
+phib = phib(:)/norm;
+psib = psib(:)/norm;
+figure;
+subplot(1,2,1);
+plot(b_angle/pi,real(phib),'b'); hold on;
+plot(b_angle/pi,imag(phib),'r'); hold on;
+plot(b_angle/pi, abs(phib),'k'); hold on;
+xlabel('$\chi/\pi$'); ylabel('$\phi(\chi)/|\phi(0)|$');
+ title('GBMS');
+subplot(1,2,2)
+plot(b_angle/pi,real(psib),'b'); hold on;
+plot(b_angle/pi,imag(psib),'r'); hold on;
+plot(b_angle/pi, abs(psib),'k'); hold on;
+xlabel('$\chi/\pi$'); ylabel('$\psi(\chi)/|\phi(0)|$');
+end
diff --git a/matlab/plot/show_geometry.m b/matlab/plot/show_geometry.m
index 5acdb7a05a9e68468698880b62b76e8f01667f84..dca0292767ff0bc27b4e710544b2db1222c6d708 100644
--- a/matlab/plot/show_geometry.m
+++ b/matlab/plot/show_geometry.m
@@ -31,7 +31,7 @@ Xfl = @(z) (R+a*cos(z)).*cos(q*z);
Yfl = @(z) (R+a*cos(z)).*sin(q*z);
Zfl = @(z) a*sin(z);
Rvec= @(z) [Xfl(z); Yfl(z); Zfl(z)];
-% xvec
+% xvec shearless
xX = @(z) (Xfl(z)-R*cos(q*z))./sqrt((Xfl(z)-R*cos(q*z)).^2+(Yfl(z)-R*sin(q*z)).^2+Zfl(z).^2);
xY = @(z) (Yfl(z)-R*sin(q*z))./sqrt((Xfl(z)-R*cos(q*z)).^2+(Yfl(z)-R*sin(q*z)).^2+Zfl(z).^2);
xZ = @(z) Zfl(z)./sqrt((Xfl(z)-R*cos(q*z)).^2+(Yfl(z)-R*sin(q*z)).^2+Zfl(z).^2);
@@ -65,12 +65,26 @@ for it_ = 1:numel(OPTIONS.TIME)
subplot(1,numel(OPTIONS.TIME),it_)
%plot magnetic geometry
if OPTIONS.PLT_MTOPO
- magnetic_topo=surf(x_tor, y_tor, z_tor); hold on;alpha 1.0;%light('Position',[-1 1 1],'Style','local')
- set(magnetic_topo,'edgecolor',[1 1 1]*0.7,'facecolor','none')
+ magnetic_topo=surf(x_tor, y_tor, z_tor); hold on;alpha 0.5;%light('Position',[-1 1 1],'Style','local')
+ set(magnetic_topo,'edgecolor',[1 1 1]*0.8,'facecolor','none')
+% set(magnetic_topo,'edgecolor','none','facecolor','white')
end
%plot field line
theta = linspace(-Nturns*pi, Nturns*pi, 512) ; % Poloidal angle
plot3(Xfl(theta),Yfl(theta),Zfl(theta)); hold on;
+ %plot fluxe tube
+ if OPTIONS.PLT_FTUBE
+ theta = linspace(-Nturns*pi, Nturns*pi, 64) ; % Poloidal angle
+ %store the shifts in an order (top left to bottom right)
+ s_x = r_o_R*[DATA.x(1) DATA.x(end) DATA.x(1) DATA.x(end)];
+ s_y = r_o_R*[DATA.y(1) DATA.y(1) DATA.y(end) DATA.y(end)];
+ for i_ = 1:4
+ vx_ = Xfl(theta) + s_x(i_)*xX(theta) + s_y(i_)*yX(theta);
+ vy_ = Yfl(theta) + s_x(i_)*xY(theta) + s_y(i_)*yY(theta);
+ vz_ = Zfl(theta) + s_x(i_)*xZ(theta) + s_y(i_)*yZ(theta);
+ plot3(vx_,vy_,vz_,'-','color',[1.0 0.6 0.6]*0.8,'linewidth',1.5); hold on;
+ end
+ end
%plot vector basis
theta = DATA.z ; % Poloidal angle
plot3(Xfl(theta),Yfl(theta),Zfl(theta),'ok'); hold on;
@@ -105,6 +119,6 @@ end
%%
axis equal
view([1,-2,1])
-
+ grid on
end
diff --git a/matlab/setup.m b/matlab/setup.m
index f5062e5787b1ab8ea2665e5bc5de8fee4aeb80fc..cc4959d6f72e1656a4916fbc91513a7e6f66a578 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -43,10 +43,10 @@ MODEL.q_e =-1.0;
MODEL.q_i = 1.0;
if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end;
% gradients L_perp/L_x
-MODEL.K_N = K_N;
-MODEL.ETA_N = ETA_N;
-MODEL.K_T = K_T;
-MODEL.ETA_T = ETA_T;
+MODEL.K_Ni = K_Ni;
+MODEL.K_Ne = K_Ne;
+MODEL.K_Ti = K_Ti;
+MODEL.K_Te = K_Te;
MODEL.GradB = GRADB; % Magnetic gradient
MODEL.CurvB = CURVB; % Magnetic curvature
MODEL.lambdaD = LAMBDAD;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 291fd2e6dfc2ca23714423e81236b01bed823c9c..69fe008dc0338f5d324827ab59ed1dbbde32874e 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -69,10 +69,10 @@ fprintf(fid,[' sigma_e = ', num2str(MODEL.sigma_e),'\n']);
fprintf(fid,[' sigma_i = ', num2str(MODEL.sigma_i),'\n']);
fprintf(fid,[' q_e = ', num2str(MODEL.q_e),'\n']);
fprintf(fid,[' q_i = ', num2str(MODEL.q_i),'\n']);
-fprintf(fid,[' K_N = ', num2str(MODEL.K_N),'\n']);
-fprintf(fid,[' ETA_N = ', num2str(MODEL.ETA_N),'\n']);
-fprintf(fid,[' K_T = ', num2str(MODEL.K_T),'\n']);
-fprintf(fid,[' ETA_T = ', num2str(MODEL.ETA_T),'\n']);
+fprintf(fid,[' K_Ne = ', num2str(MODEL.K_Ne),'\n']);
+fprintf(fid,[' K_Ni = ', num2str(MODEL.K_Ni),'\n']);
+fprintf(fid,[' K_Te = ', num2str(MODEL.K_Te),'\n']);
+fprintf(fid,[' K_Ti = ', num2str(MODEL.K_Ti),'\n']);
fprintf(fid,[' GradB = ', num2str(MODEL.GradB),'\n']);
fprintf(fid,[' CurvB = ', num2str(MODEL.CurvB),'\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
diff --git a/wk/KBM.m b/wk/KBM.m
deleted file mode 100644
index f5b5215d9c349266b6cb737412e639e3e62a6b49..0000000000000000000000000000000000000000
--- a/wk/KBM.m
+++ /dev/null
@@ -1,188 +0,0 @@
-%% QUICK RUN SCRIPT
-% This script create a directory in /results and run a simulation directly
-% from matlab framework. It is meant to run only small problems in linear
-% for benchmark and debugging purpose since it makes matlab "busy"
-%
-% SIMID = 'test_circular_geom'; % Name of the simulation
-% SIMID = 'linear_CBC'; % Name of the simulation
-SIMID = 'KBM'; % Name of the simulation
-RUN = 1; % To run or just to load
-addpath(genpath('../matlab')) % ... add
-default_plots_options
-HELAZDIR = '/home/ahoffman/HeLaZ/';
-EXECNAME = 'helaz3';
-% EXECNAME = 'helaz3_dbg';
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Set Up parameters
-CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% PHYSICAL PARAMETERS
-NU = 0.05; % Collision frequency
-TAU = 1.0; % e/i temperature ratio
-K_N = 3.0; % ele Density gradient drive
-ETA_N = 3.0/K_N; % ion Density gradient drive
-K_T = 8.0; % Temperature '''
-ETA_T = 4.5/K_T; % Temperature '''
-SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-% SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
-BETA = 0.03; % electron plasma beta
-%% GRID PARAMETERS
-P = 4;
-J = P/2;
-PMAXE = P; % Hermite basis size of electrons
-JMAXE = J; % Laguerre "
-PMAXI = P; % " ions
-JMAXI = J; % "
-NX = 16; % real space x-gridpoints
-NY = 2; % '' y-gridpoints
-LX = 2*pi/0.1; % Size of the squared frequency domain
-LY = 2*pi/0.25; % Size of the squared frequency domain
-NZ = 16; % number of perpendicular planes (parallel grid)
-NPOL = 1;
-SG = 0; % Staggered z grids option
-%% GEOMETRY
-% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
-GEOMETRY= 's-alpha';
-% GEOMETRY= 'circular';
-Q0 = 1.4; % safety factor
-SHEAR = 0.8; % magnetic shear
-NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
-EPS = 0.18; % inverse aspect ratio
-%% TIME PARMETERS
-TMAX = 25; % Maximal time unit
-DT = 5e-3; % Time step
-SPS0D = 1; % Sampling per time unit for 2D arrays
-SPS2D = 0; % Sampling per time unit for 2D arrays
-SPS3D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 1/5; % Sampling per time unit for 5D arrays
-SPSCP = 0; % Sampling per time unit for checkpoints
-JOB2LOAD= -1;
-%% OPTIONS
-LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
-% Collision operator
-% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO = 'DG';
-GKCO = 0; % gyrokinetic operator
-ABCO = 1; % interspecies collisions
-INIT_ZF = 0; ZF_AMP = 0.0;
-CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
-NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
-KERN = 0; % Kernel model (0 : GK)
-INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
-%% OUTPUTS
-W_DOUBLE = 1;
-W_GAMMA = 1; W_HF = 1;
-W_PHI = 1; W_NA00 = 1;
-W_DENS = 1; W_TEMP = 1;
-W_NAPJ = 1; W_SAPJ = 0;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% unused
-HD_CO = 0.0; % Hyper diffusivity cutoff ratio
-MU = 0.0; % Hyperdiffusivity coefficient
-INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
-MU_X = MU; %
-MU_Y = MU; %
-N_HD = 4;
-MU_Z = 1.0; %
-MU_P = 0.0; %
-MU_J = 0.0; %
-LAMBDAD = 0.0;
-NOISE0 = 1.0e-5; % Init noise amplitude
-BCKGD0 = 0.0; % Init background
-GRADB = 1.0;
-CURVB = 1.0;
-%%-------------------------------------------------------------------------
-%% RUN
-setup
-% system(['rm fort*.90']);
-% Run linear simulation
-if RUN
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
-end
-
-%% Load results
-%%
-filename = [SIMID,'/',PARAMS,'/'];
-LOCALDIR = [HELAZDIR,'results/',filename,'/'];
-% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 00;
-data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
-
-%% Short analysis
-if 1
-%% linear growth rate (adapted for 2D zpinch and fluxtube)
-trange = [0.5 1]*data.Ts3D(end);
-nplots = 3;
-lg = compute_fluxtube_growth_rate(data,trange,nplots);
-[gmax, kmax] = max(lg.g_ky(:,end));
-[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
-msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
-msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
-end
-
-if 1
-%% Ballooning plot
-options.time_2_plot = [120];
-options.kymodes = [0.1];
-options.normalized = 1;
-% options.field = 'phi';
-fig = plot_ballooning(data,options);
-end
-
-if 0
-%% Hermite-Laguerre spectrum
-% options.TIME = 'avg';
-options.P2J = 1;
-options.ST = 1;
-options.PLOT_TYPE = 'space-time';
-% options.PLOT_TYPE = 'Tavg-1D';
-% options.PLOT_TYPE = 'Tavg-2D';
-options.NORMALIZED = 0;
-options.JOBNUM = 0;
-options.TIME = [0 50];
-options.specie = 'i';
-options.compz = 'avg';
-fig = show_moments_spectrum(data,options);
-% fig = show_napjz(data,options);
-save_figure(data,fig)
-end
-
-if 0
-%% linear growth rate for 3D Zpinch (kz fourier transform)
-trange = [0.5 1]*data.Ts3D(end);
-options.keq0 = 1; % chose to plot planes at k=0 or max
-options.kxky = 1;
-options.kzkx = 0;
-options.kzky = 0;
-[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
-save_figure(data,fig)
-end
-if 0
-%% Mode evolution
-options.NORMALIZED = 0;
-options.K2PLOT = 1;
-options.TIME = [0:1000];
-options.NMA = 1;
-options.NMODES = 1;
-options.iz = 'avg';
-fig = mode_growth_meter(data,options);
-save_figure(data,fig,'.png')
-end
-
-
-if 0
-%% RH TEST
-ikx = 2; t0 = 0; t1 = data.Ts3D(end);
-[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));
-plt = @(x) squeeze(mean(real(x(1,ikx,:,it0:it1)),3))./squeeze(mean(real(x(1,ikx,:,it0)),3));
-figure
-plot(data.Ts3D(it0:it1), plt(data.PHI));
-xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
-title(sprintf('$k_x=$%2.2f, $k_y=0.00$',data.kx(ikx)))
-end
diff --git a/wk/RH_test.m b/wk/RH_test.m
deleted file mode 100644
index 463a73af7eb9d901f65dfa9481655fa90f817e2d..0000000000000000000000000000000000000000
--- a/wk/RH_test.m
+++ /dev/null
@@ -1,188 +0,0 @@
-%% QUICK RUN SCRIPT
-% This script create a directory in /results and run a simulation directly
-% from matlab framework. It is meant to run only small problems in linear
-% for benchmark and debugging purpose since it makes matlab "busy"
-%
-% SIMID = 'test_circular_geom'; % Name of the simulation
-% SIMID = 'linear_CBC'; % Name of the simulation
-SIMID = 'RH_test'; % Name of the simulation
-RUN = 1; % To run or just to load
-addpath(genpath('../matlab')) % ... add
-default_plots_options
-HELAZDIR = '/home/ahoffman/HeLaZ/';
-EXECNAME = 'helaz3';
-% EXECNAME = 'helaz3_dbg';
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Set Up parameters
-CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% PHYSICAL PARAMETERS
-NU = 0.0; % Collision frequency
-TAU = 1.0; % e/i temperature ratio
-K_N = 0.0; % ele Density gradient drive
-ETA_N = 0.0/K_N; % ion Density gradient drive
-K_T = 0.0; % Temperature '''
-ETA_T = 0.0/K_T; % Temperature '''
-SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-% SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
-BETA = 0.0; % electron plasma beta
-%% GRID PARAMETERS
-P = 20;
-J = P/2;
-PMAXE = P; % Hermite basis size of electrons
-JMAXE = J; % Laguerre "
-PMAXI = P; % " ions
-JMAXI = J; % "
-NX = 2; % real space x-gridpoints
-NY = 2; % '' y-gridpoints
-LX = 2*pi/0.05; % Size of the squared frequency domain
-LY = 2*pi/0.25; % Size of the squared frequency domain
-NZ = 16; % number of perpendicular planes (parallel grid)
-NPOL = 1;
-SG = 0; % Staggered z grids option
-%% GEOMETRY
-% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
-GEOMETRY= 's-alpha';
-% GEOMETRY= 'circular';
-Q0 = 1.4; % safety factor
-SHEAR = 0.0; % magnetic shear
-NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
-EPS = 0.18; % inverse aspect ratio
-%% TIME PARMETERS
-TMAX = 100; % Maximal time unit
-DT = 5e-3; % Time step
-SPS0D = 1; % Sampling per time unit for 2D arrays
-SPS2D = 0; % Sampling per time unit for 2D arrays
-SPS3D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 1/5; % Sampling per time unit for 5D arrays
-SPSCP = 0; % Sampling per time unit for checkpoints
-JOB2LOAD= -1;
-%% OPTIONS
-LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
-% Collision operator
-% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO = 'DG';
-GKCO = 0; % gyrokinetic operator
-ABCO = 1; % interspecies collisions
-INIT_ZF = 0; ZF_AMP = 0.0;
-CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
-NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
-KERN = 0; % Kernel model (0 : GK)
-INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom
-%% OUTPUTS
-W_DOUBLE = 1;
-W_GAMMA = 1; W_HF = 1;
-W_PHI = 1; W_NA00 = 1;
-W_DENS = 1; W_TEMP = 1;
-W_NAPJ = 1; W_SAPJ = 0;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% unused
-HD_CO = 0.0; % Hyper diffusivity cutoff ratio
-MU = 0.0; % Hyperdiffusivity coefficient
-INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
-MU_X = MU; %
-MU_Y = MU; %
-N_HD = 4;
-MU_Z = 1.0; %
-MU_P = 0.0; %
-MU_J = 0.0; %
-LAMBDAD = 0.0;
-NOISE0 = 1.0e-5; % Init noise amplitude
-BCKGD0 = 1.0; % Init background
-GRADB = 1.0;
-CURVB = 1.0;
-%%-------------------------------------------------------------------------
-%% RUN
-setup
-% system(['rm fort*.90']);
-% Run linear simulation
-if RUN
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
-end
-
-%% Load results
-%%
-filename = [SIMID,'/',PARAMS,'/'];
-LOCALDIR = [HELAZDIR,'results/',filename,'/'];
-% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 00;
-data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
-
-%% Short analysis
-if 0
-%% linear growth rate (adapted for 2D zpinch and fluxtube)
-trange = [0.5 1]*data.Ts3D(end);
-nplots = 3;
-lg = compute_fluxtube_growth_rate(data,trange,nplots);
-[gmax, kmax] = max(lg.g_ky(:,end));
-[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
-msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
-msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
-end
-
-if 0
-%% Ballooning plot
-options.time_2_plot = [120];
-options.kymodes = [0.1];
-options.normalized = 1;
-% options.field = 'phi';
-fig = plot_ballooning(data,options);
-end
-
-if 0
-%% Hermite-Laguerre spectrum
-% options.TIME = 'avg';
-options.P2J = 1;
-options.ST = 1;
-options.PLOT_TYPE = 'space-time';
-% options.PLOT_TYPE = 'Tavg-1D';
-% options.PLOT_TYPE = 'Tavg-2D';
-options.NORMALIZED = 0;
-options.JOBNUM = 0;
-options.TIME = [0 50];
-options.specie = 'i';
-options.compz = 'avg';
-fig = show_moments_spectrum(data,options);
-% fig = show_napjz(data,options);
-save_figure(data,fig)
-end
-
-if 0
-%% linear growth rate for 3D Zpinch (kz fourier transform)
-trange = [0.5 1]*data.Ts3D(end);
-options.keq0 = 1; % chose to plot planes at k=0 or max
-options.kxky = 1;
-options.kzkx = 0;
-options.kzky = 0;
-[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
-save_figure(data,fig)
-end
-if 0
-%% Mode evolution
-options.NORMALIZED = 0;
-options.K2PLOT = 1;
-options.TIME = [0:1000];
-options.NMA = 1;
-options.NMODES = 1;
-options.iz = 'avg';
-fig = mode_growth_meter(data,options);
-save_figure(data,fig,'.png')
-end
-
-
-if 1
-%% RH TEST
-ikx = 2; t0 = 0; t1 = data.Ts3D(end);
-[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));
-plt = @(x) squeeze(mean(real(x(1,ikx,:,it0:it1)),3))./squeeze(mean(real(x(1,ikx,:,it0)),3));
-figure
-plot(data.Ts3D(it0:it1), plt(data.PHI));
-xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
-title(sprintf('$k_x=$%2.2f, $k_y=0.00$',data.kx(ikx)))
-end
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index fbfbd936d360f9a9174f02915d447e2389cfe3ff..cc348d028b9958709ef7b9edc0c7880f048f5c76 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -9,7 +9,7 @@ MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
system(['mkdir -p ',MISCDIR]);
system(['mkdir -p ',LOCALDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
-% system(CMD);
+system(CMD);
% Load outputs from jobnummin up to jobnummax
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
data.localdir = LOCALDIR;
@@ -23,8 +23,8 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
% data.scale = 1;%/(data.Nx*data.Ny)^2;
-options.TAVG_0 = 250;%0.4*data.Ts3D(end);
-options.TAVG_1 = 350;%0.9*data.Ts3D(end); % Averaging times duration
+options.TAVG_0 = 400;%0.4*data.Ts3D(end);
+options.TAVG_1 = 600;%0.9*data.Ts3D(end); % Averaging times duration
options.NCUT = 4; % Number of cuts for averaging and error estimation
options.NMVA = 1; % Moving average for time traces
% options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
@@ -45,32 +45,32 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
-options.NAME = '\phi';
-% options.NAME = 'N_i^{00}';
+% options.NAME = '\phi';
+options.NAME = 'N_i^{00}';
% options.NAME = 'v_y';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'xz';
+options.PLAN = 'xy';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
-options.TIME = data.Ts3D;
-% options.TIME = [850:0.1:1000];
+% options.TIME = data.Ts3D;
+options.TIME = [00:1:800];
data.EPS = 0.1;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
end
-if 0
+if 1
%% 2D snapshots
% Options
-options.INTERP = 0;
+options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
-% options.NAME = '\phi';
-options.NAME = '\psi';
+options.NAME = '\phi';
+% options.NAME = '\psi';
% options.NAME = 'n_e';
% options.NAME = 'N_i^{00}';
% options.NAME = 'T_i';
@@ -80,7 +80,7 @@ options.PLAN = 'kxky';
% options.NAME 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [400 440];
+options.TIME = [100 200 500];
data.a = data.EPS * 2e3;
fig = photomaton(data,options);
% save_figure(data,fig)
@@ -88,12 +88,14 @@ end
if 0
%% 3D plot on the geometry
-options.INTERP = 1;
+options.INTERP = 0;
options.NAME = '\phi';
options.PLANES = [1:1:16];
-options.TIME = [15];
-options.PLT_MTOPO = 0;
-data.rho_o_R = 2e-3; % Sound larmor radius over Machine size ratio
+options.TIME = [30];
+options.PLT_MTOPO = 1;
+options.PLT_FTUBE = 1;
+data.EPS = 0.4;
+data.rho_o_R = 3e-3; % Sound larmor radius over Machine size ratio
fig = show_geometry(data,options);
save_figure(data,fig,'.png')
end
@@ -169,7 +171,7 @@ if 0
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 1;
-options.TIME = [0:160];
+options.TIME = [00:800];
options.NMA = 1;
options.NMODES = 1;
options.iz = 'avg';
diff --git a/wk/analysis_gbms.m b/wk/analysis_gbms.m
index e492b42078daed1b40fe49a707994c50deee5f34..d1957270e9db74c74b55bfac8da0cce64357a9b6 100644
--- a/wk/analysis_gbms.m
+++ b/wk/analysis_gbms.m
@@ -5,9 +5,15 @@
%%
% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/shearless_linear_cyclone/';
-resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/RH_test/';
+% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/new_RH_test/';
+% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/RH_test_kine/';
+resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/KBM/';
+% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/TEM/';
+% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/ITG/';
% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/linear_cyclone/';
% resdir = '/home/ahoffman/molix/';
+
+% system(['cd ',resdir,';','./gbms < parameters.in; cd /home/ahoffman/HeLaZ/wk']);
outfile = [resdir,'field.dat.h5'];
gbms_dat.Ts3D = h5read(outfile,'/data/var2d/time');
@@ -22,13 +28,20 @@ gbms_dat.Nz = numel(gbms_dat.z);
dky = min(gbms_dat.ky(gbms_dat.ky>0)); Ly =0;% 2*pi/dky;
gbms_dat.y = linspace(-Ly/2,Ly/2,gbms_dat.Ny+1); gbms_dat.y = gbms_dat.y(1:end-1);
gbms_dat.x = 0;
+gbms_dat.BETA = h5readatt(outfile,'/data/input','betae');
+gbms_dat.SHEAR = h5readatt(outfile,'/data/input','magnetic shear');
gbms_dat.PHI = zeros(gbms_dat.Ny,gbms_dat.Nx,gbms_dat.Nz,gbms_dat.Nt);
+gbms_dat.PSI = zeros(gbms_dat.Ny,gbms_dat.Nx,gbms_dat.Nz,gbms_dat.Nt);
gbms_dat.param_title = 'GBMS';
for it = 1:gbms_dat.Nt
tmp = h5read(outfile,['/data/var2d/phi/',sprintf('%.6d',it-1)]);
gbms_dat.PHI(:,:,:,it) = permute(tmp.real + 1i * tmp.imaginary,[2 1 3]);
-
+ if gbms_dat.BETA > 0
+ tmp = h5read(outfile,['/data/var2d/psi/',sprintf('%.6d',it-1)]);
+ gbms_dat.PSI(:,:,:,it) = permute(tmp.real + 1i * tmp.imaginary,[2 1 3]);
+ end
+
end
gbms_dat.localdir = resdir;
@@ -66,27 +79,36 @@ end
if 0
%% linear growth rate for 3D fluxtube
-trange = [10 200];
-nplots = 1;
+trange = [0.5 1]*gbms_dat.Ts3D(end);
+nplots = 3;
lg = compute_fluxtube_growth_rate(gbms_dat,trange,nplots);
+[gmax, kmax] = max(lg.g_ky(:,end));
+[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
+msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
+msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
end
-if 0
+if 1
%% Ballooning plot
-options.time_2_plot = data.Ts3D(end);
-options.kymodes = [0.5];
-options.normalized = 1;
-options.sheared = 0;
-options.field = 'phi';
-fig = plot_ballooning(gbms_dat,options);
+% options.time_2_plot = data.Ts3D(end);
+% options.kymodes = [0.5];
+% options.normalized = 1;
+% options.sheared = 0;
+% options.field = 'phi';
+% fig = plot_ballooning(gbms_dat,options);
+plot_gbms_ballooning(outfile);
+
+% plot(b_angle,phib_real); hold on;
+
end
-if 1
+if 0
%% RH TEST
-ikx = 1;
-plt = @(x) squeeze(mean(real(x(1,ikx,:,:)),3))./squeeze(mean(real(x(1,ikx,:,1)),3));
+ikx = 1; iky = 1; t0 = 0; t1 = gbms_dat.Ts3D(end);
+[~, it0] = min(abs(t0-gbms_dat.Ts3D));[~, it1] = min(abs(t1-gbms_dat.Ts3D));
+plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3));%./squeeze(mean(real(x(iky,ikx,:,it0)),3));
figure
-plot(gbms_dat.Ts3D, plt(gbms_dat.PHI));
+plot(gbms_dat.Ts3D(it0:it1), plt(gbms_dat.PHI),'k');
xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
-title(sprintf('$k_x=$%2.2f, $k_y=0.00$',gbms_dat.kx(ikx)))
+title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',gbms_dat.kx(ikx),gbms_dat.ky(iky)))
end
\ No newline at end of file
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 0ae250f5c4137f19cdb3013675b175fa99850868..22459a5e3314e67f45278c73424ac438b940bc16 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -21,12 +21,13 @@ addpath(genpath([helazdir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/Z-pinch/ZP_HP_kn_1.6_HRES/';
% folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
% folder = '/misc/gene_results/CBC/128x64x16x24x12/';
-folder = '/misc/gene_results/CBC/196x96x20x32x16_01/';
-% folder = '/misc/gene_results/CBC/128x64x16x6x4/';
+% folder = '/misc/gene_results/CBC/196x96x20x32x16_02/';
+folder = '/misc/gene_results/CBC/128x64x16x6x4/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_01/';
% folder = '/misc/gene_results/CBC/KT_4.5_128x64x16x24x12_01/';
% folder = '/misc/gene_results/CBC/KT_9_128x64x16x24x12/';
% folder = '/misc/gene_results/CBC/KT_13_large_box_128x64x16x24x12/';
+% folder = '/misc/gene_results/CBC/Lapillone_Fig6/';
gene_data = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
@@ -78,7 +79,7 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'kxky';
+options.PLAN = 'xy';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 1874ce85fd72ab361cd63bc6f73a9b059939094c..a7f051d6af91614932f89399922d25c73c967371 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -41,16 +41,18 @@ helazdir = '/home/ahoffman/HeLaZ/';
% outfile = 'CBC/64x32x16x5x3';
% outfile = 'CBC/64x128x16x5x3';
% outfile = 'CBC/128x64x16x5x3';
-% outfile = 'CBC/128x96x16x3x2_Nexc_6';
+% outfile = 'CBC/96x96x16x3x2_Nexc_6';
+% outfile = 'CBC/128x96x16x3x2';
% outfile = 'CBC/192x96x16x3x2';
% outfile = 'CBC/192x96x24x13x7';
% outfile = 'CBC/kT_11_128x64x16x5x3';
% outfile = 'CBC/kT_9_256x128x16x3x2';
-% outfile = 'CBC/kT_4.5_128x64x16x13x2';
+% outfile = 'CBC/kT_4.5_128x64x16x13x3';
% outfile = 'CBC/kT_4.5_192x96x24x13x7';
+% outfile = 'CBC/kT_4.5_128x64x16x13x7';
+% outfile = 'CBC/kT_4.5_128x96x24x15x5';
% outfile = 'CBC/kT_5.3_192x96x24x13x7';
% outfile = 'CBC/kT_13_large_box_128x64x16x5x3';
-% outfile = 'CBC/kT_13_96x96x16x3x2_Nexc_6';
% outfile = 'CBC/kT_11_96x64x16x5x3_ky_0.02';
% outfile = 'CBC/kT_scan_128x64x16x5x3';