save the scripts

fort_example.90

@@ -68,8 +68,8 @@
   K_Te    = 6.96                  !electron temperature gradient intensity
   K_Ni    = 2.22                  !ion density gradient intensity
   K_Ti    = 6.96                  !ion density temperature intensity
-  GradB     = 1                   !magnetic field gradient strength
-  CurvB     = 1                   !magnetic curvature strength
+  k_gB     = 1                   !magnetic field gradient strength
+  k_cB     = 1                   !magnetic curvature strength
   lambdaD = 0                     !Debye length (not tested when non zero)
 /
 &COLLISION_PAR

matlab/compile_results.m

@@ -160,13 +160,15 @@ while(CONTINUE)
              Nepjz  = load_3D_data(filename, 'Nepjz');
              Nepjz_ = cat(4,Nepjz_,Nepjz); clear Nepjz
                 catch
-                end
+                 disp('Cannot load Nepjz');
+               end
             end
-%             try
+            try
             [Nipjz, Ts3D, ~] = load_3D_data(filename, 'Nipjz');
              Nipjz_ = cat(4,Nipjz_,Nipjz); clear Nipjz        
-%             catch
-%             end
+            catch
+                disp('Cannot load Nipjz');
+            end
         end
         if W_DENS
             if KIN_E

matlab/extract_fig_data.m

@@ -4,7 +4,7 @@
 % tw = [3000 4000];
 % tw = [4000 4500];
 % tw = [4500 5000];
-tw = [450 650];
+tw = [100 180];
 
 fig = gcf;
 axObjs = fig.Children;

matlab/load/load_3D_data.m

@@ -6,14 +6,24 @@ function [ data, time, dt ] = load_3D_data( filename, variablename )
     
     % Find array size by loading the first output
     tmp   = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+1,'%06d')]);
-    sz    = size(tmp.real);
+    try % check if it is complex or real
+        sz  = size(tmp.real);
+        cmpx = 1;
+    catch
+        sz  = size(tmp);
+        cmpx = 0;
+    end
     % add time dimension
     sz    = [sz numel(time)];
     data     = zeros(sz);
     
     for it = 1:numel(time)
         tmp         = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
-        data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+        if cmpx
+            data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+        else
+            data(:,:,:,it) = tmp;
+        end
     end
 
 end

matlab/load/load_gene_data.m

@@ -22,6 +22,8 @@ DATA.Ny  = DATA.Nky*2-1;
 DATA.z   = h5read([folder,coofile],'/coord/z');
 DATA.Nz  = numel(DATA.z);
 
+DATA.paramshort = [num2str(DATA.Nkx),'x',num2str(DATA.Nky),'x',num2str(DATA.Nz),...
+                    'x',num2str(DATA.Nvp),'x',num2str(DATA.Nmu)];
 if numel(DATA.kx)>1
     dkx = DATA.kx(2); 
 else

matlab/load/quick_gene_load.m

@@ -73,12 +73,12 @@ path = '/home/ahoffman/gene/linear_CBC_results/';
 %----------Convergence nvpar shearless CBC
 % fname = 'CBC_salpha_nz_24_nv_scan_nw_16_adiabe.txt';
 % fname = 'CBC_miller_nz_24_nv_scan_nw_16_adiabe.txt';
-fname = 'CBC_salpha_nz_24_nv_scan_nw_16_kine.txt';
+% fname = 'CBC_salpha_nz_24_nv_scan_nw_16_kine.txt';
 % fname = 'CBC_miller_nz_24_nv_scan_nw_16_kine.txt';
 %---------- CBC
 % fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
 % fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_kine.txt';
-% fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_adiabe.txt';
+fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_adiabe.txt';
 % fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
 % fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_kine.txt';
 % fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_kine.txt';

matlab/plot/plot_metric.m

-function [ fig ] = plot_metric( data, options )
+function [ fig, arrays ] = plot_metric( data, options )
 
 % names = {'Jacobian','gradxB','gradyB','gradzB','gradz_coeff',...
 %          'gxx','gxy','gxz','gyy','gyz','gzz','hatB','hatR','hatZ'};
 names = {'gxx','gxy','gxz','gyy','gyz','gzz',...
-         'hatB', 'gradxB', 'gradyB', 'gradzB',...
+         'hatB', 'dBdx', 'dBdy', 'dBdz',...
          'Jacobian','hatR','hatZ','gradz_coeff'};
 geo_arrays = zeros(2,data.Nz,numel(names));
 
@@ -40,5 +40,7 @@ if NPLOT > 0
     axis equal
     end
 end
+%outputs
+arrays = squeeze(geo_arrays(1,:,:));
 end
 

matlab/plot/plot_radial_transport_and_spacetime.m

@@ -10,8 +10,6 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE)
     Gx_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
     Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
     Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE;
-    disp(['G_x=',sprintf('%2.2e',Gx_infty_avg),'+-',sprintf('%2.2e',Gx_infty_std)]);
-    disp(['Q_x=',sprintf('%2.2e',Qx_infty_avg),'+-',sprintf('%2.2e',Qx_infty_std)]);
 %     disp(['Q_x=',sprintf('%2.2e',Qx_infty_avg),'+-',sprintf('%2.2e',Qx_infty_std)]);
     f_avg_z      = squeeze(mean(DATA.PHI(:,:,:,:),3));
     [~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3)));
@@ -28,7 +26,7 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE)
     end
     Qx_avg    = mean(Qx_ee);
     Qx_err    =  std(Qx_ee);
-%     disp(['Q_avg=',sprintf('%2.2e',Qx_avg),'+-',sprintf('%2.2e',Qx_err)]);
+    disp(['Q_avg=',sprintf('%2.2e',Qx_avg),'+-',sprintf('%2.2e',Qx_err)]);
     %% computations
 
     % Compute zonal and non zonal energies
@@ -66,14 +64,14 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
     FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; %set(gcf, 'Position',  [500, 1000, 1000, 600])
     FIGURE.ax1 = subplot(3,1,1,'parent',FIGURE.fig);
         plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'--',...
-            'color',clr_((DATA.Pmaxi-1)/2-1,:),...
+            'color',clr_((DATA.Pmaxi-1)/2,:),...
             'DisplayName',['$\Gamma_x$ ',DATA.paramshort]); hold on;
         plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'-',...
-            'color',clr_((DATA.Pmaxi-1)/2-1,:),...
+            'color',clr_((DATA.Pmaxi-1)/2,:),...
             'DisplayName',['$Q_x$ ',DATA.paramshort]); hold on;
         ylabel('Transport')  
         if(~isnan(Qx_infty_avg))
-        plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',...
+        plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_avg, '-k',...
             'DisplayName',['$Q_{avg}=',sprintf('%2.2f',Qx_avg),'\pm',sprintf('%2.2f',Qx_err),'$']); legend('show');
             ylim([0,5*abs(Qx_infty_avg)]); 
         else

matlab/plot/show_moments_spectrum.m

@@ -2,12 +2,13 @@ function [ FIGURE ] = show_moments_spectrum( DATA, OPTIONS )
 
 Pi = DATA.Pi;
 Ji = DATA.Ji;
-if (isempty(DATA.Nipjz))
+if ~(isempty(DATA.Nipjz))
     Time_ = DATA.Ts3D;
     Nipj = sum(abs(DATA.Nipjz),3);
 else
     Time_ = DATA.Ts5D;
-    Nipj = sum(sum(sum(abs(DATA.Nipj),3),4),5);
+%     Nipj = sum(sum(sum(abs(DATA.Nipj).^2,3),4),5);
+    Nipj = sum(sum(sum(conj(DATA.Nipj).*DATA.Nipj,3),4),5);
 end
 Nipj = squeeze(Nipj);
 
@@ -89,75 +90,128 @@ if ~OPTIONS.ST
     title([TITLE,' He-La elec. spectrum']);
     end
 else
-    [JJ,PP] = meshgrid(Ji,Pi);
-    P2Ji = PP + 2*JJ;
-    % form an axis of velocity ordered moments
-    p2ji = unique(P2Ji);
-    % weights to average
-    weights = zeros(size(p2ji));
-    % space time of moments amplitude wrt to p+2j degree
-    Ni_ST = zeros(numel(p2ji),numel(Time_));
-    % fill the st diagramm by averaging every p+2j=n moments together
-    for ip = 1:numel(Pi)
-        for ij = 1:numel(Ji)
-            [~,ip2j] = min(abs(p2ji-(Pi(ip)+2*Ji(ij))));
-            Ni_ST(ip2j,:) = Ni_ST(ip2j,:) + transpose(squeeze(Nipj(ip,ij,:)));
-            weights(ip2j) = weights(ip2j) + 1;
+    if OPTIONS.P2J
+        plotname = '$\langle\sum_k |N_i^{pj}|\rangle_{p+2j=const}$';
+        [JJ,PP] = meshgrid(Ji,Pi);
+        P2Ji = PP + 2*JJ;
+        % form an axis of velocity ordered moments
+        y_ = unique(P2Ji); yname = '$p+2j$';
+        % weights to average
+        weights = zeros(size(y_));
+        % space time of moments amplitude wrt to p+2j degree
+        Ni_ST = zeros(numel(y_),numel(Time_));
+        % fill the st diagramm by averaging every p+2j=n moments together
+        for ip = 1:numel(Pi)
+            for ij = 1:numel(Ji)
+                [~,ip2j] = min(abs(y_-(Pi(ip)+2*Ji(ij))));
+                Ni_ST(ip2j,:) = Ni_ST(ip2j,:) + transpose(squeeze(Nipj(ip,ij,:)));
+                weights(ip2j) = weights(ip2j) + 1;
+            end
         end
-    end
-    % doing the average
-    for ip2j = 1:numel(p2ji)
-        Ni_ST(ip2j,:) = Ni_ST(ip2j,:)/weights(ip2j);
-    end
-    if DATA.KIN_E
-    % same for electrons!!
-    [JJ,PP] = meshgrid(Je,Pe);
-    P2Je = PP + 2*JJ;
-    % form an axis of velocity ordered moments
-    p2je = unique(P2Je);
-    % weights to average
-    weights = zeros(size(p2je));
-    % space time of moments amplitude wrt to p+2j degree
-    Ne_ST = zeros(numel(p2je),numel(Time_));
-    % fill the st diagramm by averaging every p+2j=n moments together
-    for ip = 1:numel(Pe)
-        for ij = 1:numel(Je)
-            [~,ip2j] = min(abs(p2ji-(Pe(ip)+2*Je(ij))));
-            Ne_ST(ip2j,:) = Ne_ST(ip2j,:) + transpose(squeeze(Nepj(ip,ij,:)));
-            weights(ip2j) = weights(ip2j) + 1;
+        % doing the average
+        for ip2j = 1:numel(y_)
+            Ni_ST(ip2j,:) = Ni_ST(ip2j,:)/weights(ip2j);
         end
-    end
-    % doing the average
-    for ip2j = 1:numel(p2ji)
-        Ne_ST(ip2j,:) = Ne_ST(ip2j,:)/weights(ip2j);
-    end 
+        if DATA.KIN_E
+        % same for electrons!!
+        [JJ,PP] = meshgrid(Je,Pe);
+        P2Je = PP + 2*JJ;
+        % form an axis of velocity ordered moments
+        p2je = unique(P2Je);
+        % weights to average
+        weights = zeros(size(p2je));
+        % space time of moments amplitude wrt to p+2j degree
+        Ne_ST = zeros(numel(p2je),numel(Time_));
+        % fill the st diagramm by averaging every p+2j=n moments together
+        for ip = 1:numel(Pe)
+            for ij = 1:numel(Je)
+                [~,ip2j] = min(abs(y_-(Pe(ip)+2*Je(ij))));
+                Ne_ST(ip2j,:) = Ne_ST(ip2j,:) + transpose(squeeze(Nepj(ip,ij,:)));
+                weights(ip2j) = weights(ip2j) + 1;
+            end
+        end
+        % doing the average
+        for ip2j = 1:numel(y_)
+            Ne_ST(ip2j,:) = Ne_ST(ip2j,:)/weights(ip2j);
+        end 
+        end
+        ticks_labels = {};
+    else % We just order our moments w.r.t. to the convention ..
+         % (0,0) (1,0) (2,0) (0,1) (3,0) (1,1) (4,0) (2,1) (0,2) etc.
+        plotname = '$\langle\sum_k |N_i^{pj}|^2\rangle_z$';
+        Nmoments = numel(Nipj(:,:,1)); % total number of moments
+        HL_deg   = zeros(Nmoments,2);   % list of degrees, first column Hermite, second Laguerre
+        im  = 2;
+        deg = 1; % the zero degree is always here first place
+        ticks_labels = cell(10,1);
+        ticks_labels{1} = '(0,0)';
+        while(im<=Nmoments)
+        FOUND = 1;
+            while(FOUND) % As we find a pair matching the degree we retry
+            FOUND = 0;
+            for ij = 1:DATA.Jmaxi
+            for ip = 1:DATA.Pmaxi
+                if((ip-1)+2*(ij-1) == deg)
+                    % Check if the pair is already added
+                    check_ = HL_deg == [DATA.Pi(ip) DATA.Pi(ij)];
+                    check_ = sum(check_(:,1) .* check_(:,2));
+                    if ~check_ % if not add it
+                        HL_deg(im,1) = DATA.Pi(ip);
+                        HL_deg(im,2) = DATA.Pi(ij);
+                        ticks_labels{im} = ['(',num2str(DATA.Pi(ip)),',',num2str(DATA.Ji(ij)),')'];
+                        im = im + 1; FOUND = 1;
+                    end
+                end
+                end
+            end
+            end
+            % No pair found anymore, increase the degree
+            deg = deg + 1;
+        end
+        
+        % form an axis of velocity ordered moments
+        y_ = 1:Nmoments; yname = '$(P,J)$';
+        % space time of moments amplitude wrt to p+2j degree
+        Ni_ST = zeros(Nmoments,numel(Time_));
+        for i_ = 1:Nmoments
+           Ni_ST(i_,:) = Nipj(HL_deg(i_,1)+1,HL_deg(i_,2)+1,:); 
+        end
+        if DATA.KIN_E
+        % space time of moments amplitude wrt to p+2j degree
+        Ne_ST = zeros(Nmoments,numel(Time_));
+        for i_ = 1:Nmoments
+           Ne_ST = Nepj(HL_deg(i_,1)+1,HL_deg(i_,2)+1,:); 
+        end
+        end    
+         
     end
     % plots
     % scaling
-%     plt = @(x,ip2j) x./max(max(x));
     if OPTIONS.NORMALIZED
-    plt = @(x,ip2j) x(ip2j,:)./max(x(ip2j,:));
+    plt = @(x,i) x(i,:)./max(x(i,:));
     else
-    plt = @(x,ip2j) x;
+    plt = @(x,i) x;
     end
     if DATA.KIN_E
     subplot(2,1,1)
     end
-        imagesc(Time_,p2ji,plt(Ni_ST,1:numel(p2ji))); 
-        set(gca,'YDir','normal')        
-%         pclr = pcolor(XX,YY,plt(Ni_ST));
-%         set(pclr, 'edgecolor','none'); hold on;
-    xlabel('$t$'); ylabel('$p+2j$')
-    title('$\langle\sum_k |N_i^{pj}|\rangle_{p+2j=const}$')
+    
+    imagesc(Time_,y_,plt(Ni_ST,1:numel(y_))); 
+    set(gca,'YDir','normal')        
+    xlabel('$t$'); ylabel(yname);
+    if ~isempty(ticks_labels)
+        yticks(y_);
+        yticklabels(ticks_labels)
+    end
+    title(plotname)
+    
     if DATA.KIN_E
     subplot(2,1,2)
-        imagesc(Time_,p2je,plt(Ne_ST,1:numel(p2ji))); 
+        imagesc(Time_,p2je,plt(Ne_ST,1:numel(y_))); 
         set(gca,'YDir','normal')
-%         pclr = pcolor(XX,YY,plt(Ne_ST)); 
-%         set(pclr, 'edgecolor','none'); hold on;
-    xlabel('$t$'); ylabel('$p+2j$')
-    title('$\langle\sum_k |N_e^{pj}|\rangle_{p+2j=const}$')
-    suptitle(DATA.param_title);
+        xlabel('$t$'); ylabel(yname)
+        title(plotname)
+        suptitle(DATA.param_title);
     end
 
 end

matlab/setup.m

@@ -52,8 +52,8 @@ 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.k_gB   = k_gB;      % Magnetic gradient
+MODEL.k_cB   = k_cB;      % Magnetic curvature
 MODEL.lambdaD = LAMBDAD;
 % Collision parameters
 COLL.collision_model = ['''',CO,''''];

matlab/write_fort90.m

@@ -77,8 +77,8 @@ fprintf(fid,['  K_Ne    = ', num2str(MODEL.K_Ne),'\n']);
 fprintf(fid,['  K_Te    = ', num2str(MODEL.K_Te),'\n']);
 fprintf(fid,['  K_Ni    = ', num2str(MODEL.K_Ni),'\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,['  k_gB   = ', num2str(MODEL.k_gB),'\n']);
+fprintf(fid,['  k_cB   = ', num2str(MODEL.k_cB),'\n']);
 fprintf(fid,['  lambdaD = ', num2str(MODEL.lambdaD),'\n']);
 fprintf(fid,['  beta    = ', num2str(MODEL.beta),'\n']);
 fprintf(fid,'/\n');

testcases/fort.90

@@ -63,8 +63,8 @@
   K_Ti    = 6.92
   K_Ne    = 0
   K_Te    = 0
-  GradB   = 1
-  CurvB   = 1
+  k_gB   = 1
+  k_cB   = 1
   lambdaD = 0
   beta    = 0
 /

testcases/matlab_testscripts/Hallenbert.m

@@ -24,8 +24,8 @@ J       = 2;
 Q0      = 1.0;       % safety factor
 SHEAR   = 0.0;       % magnetic shear
 EPS     = 0.0;      % inverse aspect ratio
-GRADB   = 1.0;       % Magnetic  gradient
-CURVB   = 1.0;       % Magnetic  curvature
+k_gB   = 1.0;       % Magnetic  gradient
+k_cB   = 1.0;       % Magnetic  curvature
 SG      = 0;         % Staggered z grids option
 %% TIME PARAMETERS
 TMAX    = 2000;  % Maximal time unit

testcases/matlab_testscripts/cyclone_test_case.m

@@ -23,8 +23,8 @@ J       = 2;
 Q0      = 1.4;       % safety factor
 SHEAR   = 0.0;       % magnetic shear
 EPS     = 0.18;      % inverse aspect ratio
-GRADB   = 1.0;       % Magnetic  gradient
-CURVB   = 1.0;       % Magnetic  curvature
+k_gB   = 1.0;       % Magnetic  gradient
+k_cB   = 1.0;       % Magnetic  curvature
 SG      = 1;         % Staggered z grids option
 %% TIME PARAMETERS
 TMAX    = 500;  % Maximal time unit

testcases/matlab_testscripts/linear_1D_entropy_mode.m

@@ -75,8 +75,8 @@ 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;
+k_gB   = 1.0;
+k_cB   = 1.0;
 %% PARAMETER SCANS
 
 if 1

testcases/matlab_testscripts/linear_damping.m

@@ -67,8 +67,8 @@ MU_J    = 0.0;     %
 LAMBDAD = 0.0;
 NOISE0  = 0*1.0e-5; % Init noise amplitude
 BCKGD0  = 1.0;    % Init background
-GRADB   = 0.0;
-CURVB   = 0.0;
+k_gB   = 0.0;
+k_cB   = 0.0;
 %% PARAMETER SCANS
 
 if 1

testcases/miller_example_w_salpha/fort_00.90

@@ -63,8 +63,8 @@
   K_Ti    = 6.92
   K_Ne    = 0
   K_Te    = 0
-  GradB   = 1
-  CurvB   = 1
+  k_gB   = 1
+  k_cB   = 1
   lambdaD = 0
   beta    = 0
 /

testcases/miller_example_w_salpha/fort_01.90

@@ -63,8 +63,8 @@
   K_Ti    = 6.92
   K_Ne    = 0
   K_Te    = 0
-  GradB   = 1
-  CurvB   = 1
+  k_gB   = 1
+  k_cB   = 1
   lambdaD = 0
   beta    = 0
 /

testcases/zpinch_example/fort.90

@@ -62,8 +62,8 @@
   K_Te    = 0.4
   K_Ni    = 2.0
   K_Ti    = 0.4
-  GradB   = 1
-  CurvB   = 1
+  k_gB   = 1
+  k_cB   = 1
   lambdaD = 0
   beta    = 0
 /