diff --git a/matlab/LinearFit_s.m b/matlab/LinearFit_s.m
index 6712007afeb58c49b1fadc0bd92599930c49ca77..af790195cad7d1b8c99f3a678e91220b33d5bb1f 100644
--- a/matlab/LinearFit_s.m
+++ b/matlab/LinearFit_s.m
@@ -1,4 +1,4 @@
-function [gamma,fit] = LinearFit_s(time,Na00abs)
+function [gamma,t,gammaoft] = LinearFit_s(time,Na00abs)
% LinearFit computes the growth rate and frequency from the time evolution of Napj
% - adapted from MOLI (B.J. Frei)
%
@@ -29,7 +29,6 @@ function [gamma,fit] = LinearFit_s(time,Na00abs)
gamma = mean(gammaoft); % ... take the mean of gamma over the time window
% Return gamma(t) for amplitude ratio method
- fit.gammaoft = gammaoft;
- fit.t = time;
+ t = 0.5*(time(2:end)+time(1:end-1));
end % ... end function
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 2bb4b5ef8f08bc0555c77bb66bb1d49ff977492a..5fa860f7b7b6e0c35479dd6f834a0a8dcbceaa59 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -5,7 +5,8 @@ JOBNUM = JOBNUMMIN; JOBFOUND = 0;
DATA.TJOB_SE = []; % Start and end times of jobs
DATA.NU_EVOL = []; % evolution of parameter nu between jobs
DATA.CO_EVOL = []; % evolution of CO
-DATA.MU_EVOL = []; % evolution of parameter mu between jobs
+DATA.MUx_EVOL = []; % evolution of parameter mu between jobs
+DATA.MUy_EVOL = []; % evolution of parameter mu between jobs
DATA.K_N_EVOL = []; %
DATA.L_EVOL = []; %
DATA.DT_EVOL = []; %
@@ -159,7 +160,8 @@ while(CONTINUE)
DATA.TJOB_SE = [DATA.TJOB_SE Ts0D(1) Ts0D(end)];
DATA.NU_EVOL = [DATA.NU_EVOL DATA.NU DATA.NU];
DATA.CO_EVOL = [DATA.CO_EVOL DATA.CO DATA.CO];
- DATA.MU_EVOL = [DATA.MU_EVOL DATA.MU DATA.MU];
+ DATA.MUx_EVOL = [DATA.MUx_EVOL DATA.MUx DATA.MUx];
+ DATA.MUy_EVOL = [DATA.MUy_EVOL DATA.MUy DATA.MUy];
DATA.K_N_EVOL = [DATA.K_N_EVOL DATA.K_N DATA.K_N];
DATA.L_EVOL = [DATA.L_EVOL DATA.L DATA.L];
DATA.DT_EVOL = [DATA.DT_EVOL DATA.DT_SIM DATA.DT_SIM];
@@ -206,7 +208,11 @@ DATA.Nx = Nx; DATA.Ny = Ny; DATA.Nz = Nz; DATA.Nkx = Nkx; DATA.Nky = Nky;
DATA.Pmaxe = numel(Pe); DATA.Pmaxi = numel(Pi); DATA.Jmaxe = numel(Je); DATA.Jmaxi = numel(Ji);
DATA.dir = DIRECTORY;
DATA.localdir = ['..',DIRECTORY(20:end)];
-
+DATA.param_title=['$\nu_{',DATA.CONAME,'}=$', num2str(DATA.NU), ...
+ ', $\kappa_N=$',num2str(DATA.K_N),', $L=',num2str(DATA.L),'$, $N=',...
+ num2str(DATA.Nx),'$, $(P,J)=(',num2str(DATA.PMAXI),',',...
+ num2str(DATA.JMAXI),')$,',' $\mu_{hd}=$(',num2str(DATA.MUx),...
+ ',',num2str(DATA.MUy),')'];
JOBNUM = LASTJOB;
filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
diff --git a/matlab/computeFXYZ.m b/matlab/computeFXYZ.m
new file mode 100644
index 0000000000000000000000000000000000000000..8982392941b3bcc13d3820a442e4106f544baf34
--- /dev/null
+++ b/matlab/computeFXYZ.m
@@ -0,0 +1,58 @@
+function output=computeFXYZ(gene_data,varargin)
+%Function for antitrasforming a Fourier field to real space
+%
+% INPUTS:
+% gene_data -> 3D field in Fourier
+% x_local -> Fourier in x (1st dim.)
+%
+% [output]=COMPUTEFXYZ(gene_data)
+
+%%
+if nargin==1
+ x_local=1;
+else
+ x_local=varargin{1};
+end
+
+if x_local
+ [nx,nky,nz]=size(gene_data);
+ %extend to whole ky by imposing reality condition
+ ny=2*nky-1;
+
+ if ny~=1
+ %note, we need one extra point which we set to zero for the ifft
+ spectrumKxKyZ=zeros(nx,ny,nz);
+ spectrumKxKyZ(:,1:nky,:)=gene_data(:,:,:);
+ spectrumKxKyZ(1,(nky+1):(ny),:)=conj(gene_data(1,nky:-1:2,:));
+ spectrumKxKyZ(2:nx,(nky+1):(ny),:)=conj(gene_data(nx:-1:2,nky:-1:2,:));
+ else
+ %pad with zeros to interpolate on fine scale
+ ny=20;
+ spectrumKxKyZ=zeros(nx,ny,nz);
+ spectrumKxKyZ(:,2,:)=gene_data(:,:,:);
+ end
+ %inverse fft, symmetric as we are using real data
+ spectrumXKyZ=nx*ifft(spectrumKxKyZ,[],1);
+ output=ny*ifft(spectrumXKyZ,[],2,'symmetric');
+ clear spectrumKxKyZ
+
+else
+
+ [nx,nky,nz]=size(gene_data);
+ % extend to whole ky by imposingreality condition
+ ny=2*nky-1;
+ if ny~=1
+ spectrumXKyZ=zeros(nx,ny,nz);
+ spectrumXKyZ(:,1:nky,:)=gene_data(:,:,:);
+ spectrumXKyZ(1,(nky+1):(ny),:)=conj(gene_data(1,nky:-1:2,:));
+ spectrumXKyZ(2:nx,(nky+1):(ny),:)=conj(gene_data(2:1:nx,nky:-1:2,:));
+ else
+ %pad with zeros to interpolate on fine scale
+ ny=50;
+ spectrumXKyZ=zeros(nx,ny,nz);
+ spectrumXKyZ(:,2,:)=gene_data(:,:,:);
+
+ end
+ output=ny*ifft(spectrumXKyZ,[],2,'symmetric');
+ clear spectrumK=XKyZ
+end
diff --git a/matlab/default_plots_options.m b/matlab/default_plots_options.m
index f8d5496697c9d49cedc9d9c9a7dd7aa11b822cc2..4654b0bf2771412297353b63245ef0601de3b6c2 100644
--- a/matlab/default_plots_options.m
+++ b/matlab/default_plots_options.m
@@ -9,3 +9,16 @@ line_colors = [[0, 0.4470, 0.7410];[0.8500, 0.3250, 0.0980];[0.9290, 0.6940, 0.1
line_styles = {'-','-.','--',':'};
marker_styles = {'^','v','o','s'};
latexize = @(x) ['$',x,'$'];
+if 0
+%%
+set(0,'defaulttextInterpreter','tex')
+set(0, 'defaultAxesTickLabelInterpreter','tex');
+set(0, 'defaultLegendInterpreter','tex');
+set(0,'defaultAxesFontSize',16)
+set(0, 'DefaultLineLineWidth', 2.0);
+line_colors = [[0, 0.4470, 0.7410];[0.8500, 0.3250, 0.0980];[0.9290, 0.6940, 0.1250];...
+ [0.4940, 0.1840, 0.5560];[0.4660, 0.6740, 0.1880];[0.3010, 0.7450, 0.9330];[0.6350, 0.0780, 0.1840]];
+line_styles = {'-','-.','--',':'};
+marker_styles = {'^','v','o','s'};
+latexize = @(x) ['$',x,'$'];
+end
\ No newline at end of file
diff --git a/matlab/distinguishable_colors.m b/matlab/distinguishable_colors.m
new file mode 100644
index 0000000000000000000000000000000000000000..eb1e24ee922f84c221b75709e90673462dce8144
--- /dev/null
+++ b/matlab/distinguishable_colors.m
@@ -0,0 +1,147 @@
+function colors = distinguishable_colors(n_colors,bg,func)
+% DISTINGUISHABLE_COLORS: pick colors that are maximally perceptually distinct
+%
+% When plotting a set of lines, you may want to distinguish them by color.
+% By default, Matlab chooses a small set of colors and cycles among them,
+% and so if you have more than a few lines there will be confusion about
+% which line is which. To fix this problem, one would want to be able to
+% pick a much larger set of distinct colors, where the number of colors
+% equals or exceeds the number of lines you want to plot. Because our
+% ability to distinguish among colors has limits, one should choose these
+% colors to be "maximally perceptually distinguishable."
+%
+% This function generates a set of colors which are distinguishable
+% by reference to the "Lab" color space, which more closely matches
+% human color perception than RGB. Given an initial large list of possible
+% colors, it iteratively chooses the entry in the list that is farthest (in
+% Lab space) from all previously-chosen entries. While this "greedy"
+% algorithm does not yield a global maximum, it is simple and efficient.
+% Moreover, the sequence of colors is consistent no matter how many you
+% request, which facilitates the users' ability to learn the color order
+% and avoids major changes in the appearance of plots when adding or
+% removing lines.
+%
+% Syntax:
+% colors = distinguishable_colors(n_colors)
+% Specify the number of colors you want as a scalar, n_colors. This will
+% generate an n_colors-by-3 matrix, each row representing an RGB
+% color triple. If you don't precisely know how many you will need in
+% advance, there is no harm (other than execution time) in specifying
+% slightly more than you think you will need.
+%
+% colors = distinguishable_colors(n_colors,bg)
+% This syntax allows you to specify the background color, to make sure that
+% your colors are also distinguishable from the background. Default value
+% is white. bg may be specified as an RGB triple or as one of the standard
+% "ColorSpec" strings. You can even specify multiple colors:
+% bg = {'w','k'}
+% or
+% bg = [1 1 1; 0 0 0]
+% will only produce colors that are distinguishable from both white and
+% black.
+%
+% colors = distinguishable_colors(n_colors,bg,rgb2labfunc)
+% By default, distinguishable_colors uses the image processing toolbox's
+% color conversion functions makecform and applycform. Alternatively, you
+% can supply your own color conversion function.
+%
+% Example:
+% c = distinguishable_colors(25);
+% figure
+% image(reshape(c,[1 size(c)]))
+%
+% Example using the file exchange's 'colorspace':
+% func = @(x) colorspace('RGB->Lab',x);
+% c = distinguishable_colors(25,'w',func);
+% Copyright 2010-2011 by Timothy E. Holy
+ % Parse the inputs
+ if (nargin < 2)
+ bg = [1 1 1]; % default white background
+ else
+ if iscell(bg)
+ % User specified a list of colors as a cell aray
+ bgc = bg;
+ for i = 1:length(bgc)
+ bgc{i} = parsecolor(bgc{i});
+ end
+ bg = cat(1,bgc{:});
+ else
+ % User specified a numeric array of colors (n-by-3)
+ bg = parsecolor(bg);
+ end
+ end
+
+ % Generate a sizable number of RGB triples. This represents our space of
+ % possible choices. By starting in RGB space, we ensure that all of the
+ % colors can be generated by the monitor.
+ n_grid = 30; % number of grid divisions along each axis in RGB space
+ x = linspace(0,1,n_grid);
+ [R,G,B] = ndgrid(x,x,x);
+ rgb = [R(:) G(:) B(:)];
+ if (n_colors > size(rgb,1)/3)
+ error('You can''t readily distinguish that many colors');
+ end
+
+ % Convert to Lab color space, which more closely represents human
+ % perception
+ if (nargin > 2)
+ lab = func(rgb);
+ bglab = func(bg);
+ else
+ C = makecform('srgb2lab');
+ lab = applycform(rgb,C);
+ bglab = applycform(bg,C);
+ end
+ % If the user specified multiple background colors, compute distances
+ % from the candidate colors to the background colors
+ mindist2 = inf(size(rgb,1),1);
+ for i = 1:size(bglab,1)-1
+ dX = bsxfun(@minus,lab,bglab(i,:)); % displacement all colors from bg
+ dist2 = sum(dX.^2,2); % square distance
+ mindist2 = min(dist2,mindist2); % dist2 to closest previously-chosen color
+ end
+
+ % Iteratively pick the color that maximizes the distance to the nearest
+ % already-picked color
+ colors = zeros(n_colors,3);
+ lastlab = bglab(end,:); % initialize by making the "previous" color equal to background
+ for i = 1:n_colors
+ dX = bsxfun(@minus,lab,lastlab); % displacement of last from all colors on list
+ dist2 = sum(dX.^2,2); % square distance
+ mindist2 = min(dist2,mindist2); % dist2 to closest previously-chosen color
+ [~,index] = max(mindist2); % find the entry farthest from all previously-chosen colors
+ colors(i,:) = rgb(index,:); % save for output
+ lastlab = lab(index,:); % prepare for next iteration
+ end
+end
+function c = parsecolor(s)
+ if ischar(s)
+ c = colorstr2rgb(s);
+ elseif isnumeric(s) && size(s,2) == 3
+ c = s;
+ else
+ error('MATLAB:InvalidColorSpec','Color specification cannot be parsed.');
+ end
+end
+function c = colorstr2rgb(c)
+ % Convert a color string to an RGB value.
+ % This is cribbed from Matlab's whitebg function.
+ % Why don't they make this a stand-alone function?
+ rgbspec = [1 0 0;0 1 0;0 0 1;1 1 1;0 1 1;1 0 1;1 1 0;0 0 0];
+ cspec = 'rgbwcmyk';
+ k = find(cspec==c(1));
+ if isempty(k)
+ error('MATLAB:InvalidColorString','Unknown color string.');
+ end
+ if k~=3 || length(c)==1,
+ c = rgbspec(k,:);
+ elseif length(c)>2,
+ if strcmpi(c(1:3),'bla')
+ c = [0 0 0];
+ elseif strcmpi(c(1:3),'blu')
+ c = [0 0 1];
+ else
+ error('MATLAB:UnknownColorString', 'Unknown color string.');
+ end
+ end
+end
\ No newline at end of file
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index fa28a51fe520be2f73dfaa26d93940bf92844994..1dd71204e95f7d0c1a756d4778c733587667c2e7 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -4,7 +4,7 @@ dataObjs = axObjs.Children;
X_ = dataObjs(1).XData;
Y_ = dataObjs(1).YData;
-n0 = 1500;
+n0 = 0001;
figure;
plot(X_(n0:end),Y_(n0:end));
plot(X_(n0:end)-X_(n0),Y_(n0:end));
\ No newline at end of file
diff --git a/matlab/ifourier_GENE.m b/matlab/ifourier_GENE.m
new file mode 100644
index 0000000000000000000000000000000000000000..15a94dabb58267d2487838b668088f3666c12ead
--- /dev/null
+++ b/matlab/ifourier_GENE.m
@@ -0,0 +1,26 @@
+function [ field_r ] = ifourier_GENE( field_c, dims )
+%IFOURIER_FIELD_GENE method of ifourier used in GENE post processing
+% This fourier transform back from the half complex plane to a full real
+% space, in 2 or 3D (dim). Compied from computeFXYZ of GENE for
+% comparison purpose.
+sz = size(field_c);
+nkx = sz(1);
+field_r = zeros(dims);
+
+if numel(dims) == 2 %2D
+ nx = dims(1); ny = dims(2);
+ %note, we need one extra point which we set to zero for the ifft
+ spectrumKxKyZ = zeros(nx,ny);
+ spectrumKxKyZ(1:nkx,:) = field_c;
+ spectrumKxKyZ((nkx):(nx),1) = conj(field_c(nkx:-1:2));
+ spectrumKxKyZ((nkx):(nx),2:ny) = conj(field_c(nkx:-1:2,ny:-1:2));
+
+% spectrumKxYZ = ny*ifft(spectrumKxKyZ,[],2);
+% field_r = nx*ifft(spectrumKxYZ,[],1,'symmetric');
+
+ spectrumKxYZ = ifft(spectrumKxKyZ,[],2);
+ field_r = ifft(spectrumKxYZ,[],1,'symmetric');
+end
+
+end
+
diff --git a/matlab/loadGeneral_GENE.m b/matlab/loadGeneral_GENE.m
new file mode 100644
index 0000000000000000000000000000000000000000..d771c033429e318e03287ed78028901d2f4712d7
--- /dev/null
+++ b/matlab/loadGeneral_GENE.m
@@ -0,0 +1,66 @@
+function [general]= loadGeneral_GENE(folder,fileNumber,varargin)
+% Loads all the input params of a simulation and returns them in a
+% structured way, If [show] is provided then results are shown on screen
+% in a compact version.
+% To see the way the structure is filled, read the code or try it once.
+%
+% INPUTS:
+% folder -> folder cointaining the simulation data
+% fileNumber -> actual run according to our numeration
+% [show] -> boolean to disply information on screen
+% [use_h5] -> define whether h5 files are to be used (default 1)
+%
+% [general]= LOADGENERAL(folder,fileNumber,[show])
+%
+
+%%
+%we here assume that all the runs are the same in terms of parameters.
+fileNumber=fileNumber(1);
+
+use_h5=0;
+
+
+[allParams]=loadParams(folder,fileNumber,use_h5);
+
+[coord]=loadCoord(folder,fileNumber,allParams.info.x_local,use_h5,allParams.specs.s1.name);
+
+[geometry]=loadGeometry(folder,fileNumber,allParams.info.geom,use_h5,allParams.info.x_local,allParams.box.nx);
+
+
+% if nargin==3 && varargin{1}==1
+% describe_simulation(folder,fileNumber);
+% end
+
+general=struct('coord',coord, 'geometry',geometry, 'folder', folder, 'fileNumber',fileNumber, 'allParams',allParams);
+
+if ~general.allParams.info.x_local
+ for i_s=1:general.allParams.specs.ns
+ act_spec=general.allParams.specs.(genvarname(['s',num2str(i_s)])).name;
+ if use_h5
+ file=fileNameH5(folder,['profiles_',act_spec],fileNumber );
+ T=h5read(file,'/temp/T');
+ omt=h5read(file,'/temp/omt');
+ n=h5read(file,'/density/n');
+ omn=h5read(file,'/density/omn');
+ x_o_a=h5read(file,'/position/x_o_a');
+ x_o_rho_ref=h5read(file,'/position/x_o_rho_ref');
+ general.allParams.specs.(genvarname(['s',num2str(i_s)])).profiles=struct('T',T,'n',n,'omn',omn,'omt',omt,'x_o_a',x_o_a,'x_o_rho_ref',x_o_rho_ref);
+ else
+ file=fileName_std(folder,['profiles_',act_spec],fileNumber );
+ fid=fopen(file);
+ nel=numel(regexp(fgetl(fid),'\s+','split'));
+ fgetl(fid);
+ data=textscan(fid,'%f');
+ data=reshape(data{:}, [nel-1 numel(data{:})/(nel-1)]);
+ fclose(fid);
+ general.allParams.specs.(genvarname(['s',num2str(i_s)])).profiles=...
+ struct('T',data(3,:)','n',data(4,:)','omt',data(5,:)','omn',data(6,:)','x_o_a',data(1,:)','x_o_rho_ref',data(2,:)');
+
+ end
+ end
+
+ general.prefactors=compute_prefactors(general.allParams,general.geometry,numel(coord.ky),numel(coord.z));
+end
+
+
+end
diff --git a/matlab/load_field_gene.m b/matlab/load_field_gene.m
new file mode 100644
index 0000000000000000000000000000000000000000..a323a9a277240a7dc0a28dd5f6af36eff4b0256f
--- /dev/null
+++ b/matlab/load_field_gene.m
@@ -0,0 +1,98 @@
+% folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_mu_1e-2_KHR_2/';
+% folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-3_gyroLES/';
+% folder = '/misc/gene_results/HP_fig_2c_gyroLES/';
+folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_mu_1e-2_SGDK/';
+
+%% General info (domain etc.)
+varargin{1} = 0; varargin{2} = 0;
+general = loadGeneral_GENE(folder,-1,varargin);
+
+nx = general.allParams.box.nx; Lx = general.allParams.box.Lx;
+ny = general.allParams.box.ny; Ly = general.allParams.box.Ly;
+nz = general.allParams.box.nz;
+kx = general.coord.kx;
+ky = general.coord.ky;
+[KY,KX] = meshgrid(ky,kx);
+x = Lx/nx.*(-nx/2:nx/2-1);
+y = Ly/ny.*(-ny+1:ny-1);
+
+%% Load dens_i
+steps_range = 0:1:500; nt = numel(steps_range);
+time_dens = 1:nt;
+dens_xt = zeros(nt,nx);
+% figure
+for i = 1:nt
+ field_step = steps_range(i);
+
+ [field_c, t] = loadMomentumCpx(folder,-1,field_step-1,'ions','dens','none',0,[nx ny nz],0,general.allParams.data_format);
+
+ field_r = computeFXYZ(field_c,1);
+ dens_xt(i,:) = mean(field_r(:,:,1),2);
+ time_dens(i) = t;
+end
+
+if 1
+%% Load phi
+steps_range = 0:2:1000; nt = numel(steps_range);
+time_phi = 1:nt;
+phi_xt = zeros(nt,nx);
+% figure
+for i = 1:nt
+ field_step = steps_range(i);
+
+ [field_c, t] = loadFieldPhiCpx(folder,-1,field_step-1,'',0,0,[nx ny nz],general.allParams.data_format);
+
+ field_r = computeFXYZ(field_c,1);
+ phi_xt(i,:) = mean(field_r(:,:,1),2);
+ time_phi(i) = t;
+end
+end
+
+if 0
+%% Load v_ExB
+steps_range = 0:2:1000; nt = numel(steps_range);
+time_vExB = 1:nt;
+vExB_xt = zeros(nt,nx);
+% figure
+for i = 1:nt
+ field_step = steps_range(i);
+
+ [field_c, t] = loadFieldPhiCpx(folder,-1,field_step-1,'',0,0,[nx ny nz],general.allParams.data_format);
+
+ field_r = computeFXYZ(-KX.*field_c,1);
+ vExB_xt(i,:) = mean(field_r(:,:,1),2);
+ time_vExB(i) = t;
+end
+end
+if 1
+%% compute Gx
+steps_n = 0:1:500; nt1 = numel(steps_n);
+steps_v = 0:2:1000; nt2 = numel(steps_v);
+time_vExB = 1:nt1;
+Gx = zeros(nt,1);
+% figure
+for i = 1:nt1
+ field_step = steps_n(i);
+ [field_c, t] = loadMomentumCpx(folder,-1,field_step-1,'ions','dens','none',0,[nx ny nz],0,general.allParams.data_format);
+ dens = computeFXYZ(field_c,1);
+ time_dens(i) = t;
+
+ field_step = steps_v(i);
+ [field_c, t] = loadFieldPhiCpx(folder,-1,field_step-1,'',0,0,[nx ny nz],general.allParams.data_format);
+ vExB = computeFXYZ(-KY.*field_c,1);
+ time_vExB(i) = t;
+
+ Gx(i) = mean(mean(squeeze(dens(:,:,1).*vExB(:,:,1))));
+end
+end
+%% Plot
+figure;
+[X_TX,Y_TX] = meshgrid(time_dens,x);
+% pclr=pcolor(X_TX',Y_TX',dens_xt);
+pclr=pcolor(X_TX',Y_TX',phi_xt);
+set(pclr, 'edgecolor','none');
+colormap(bluewhitered(256))
+
+%%
+% figure
+% plot(time_dens,Gx)
diff --git a/matlab/load_params.m b/matlab/load_params.m
index 1d44a63a5dac5bd3c4457f27dc7337dc110c6ac5..78020b7f634d44b6650454258210d9e2fb06cfe5 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -20,6 +20,8 @@ DATA.L = h5readatt(filename,'/data/input','Lx');
DATA.CLOS = h5readatt(filename,'/data/input','CLOS');
DATA.DT_SIM = h5readatt(filename,'/data/input','dt');
DATA.MU = h5readatt(filename,'/data/input','mu');
+DATA.MUx = -99;%h5readatt(filename,'/data/input','mu_x');
+DATA.MUy = -99;%h5readatt(filename,'/data/input','mu_y');
% MU = str2num(filename(end-18:end-14)); %bad...
DATA.W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y';
DATA.W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y';
diff --git a/matlab/photomaton.m b/matlab/photomaton.m
index ef43d8740afb61bafdce67a37b4e654f9f2153dc..66b6f04832b2f43315fb7fe50570158fc60eb802 100644
--- a/matlab/photomaton.m
+++ b/matlab/photomaton.m
@@ -4,22 +4,27 @@ function [ FIGURE ] = photomaton( DATA,OPTIONS )
toplot = process_field(DATA,OPTIONS);
FNAME = toplot.FILENAME;
FRAMES = toplot.FRAMES;
-
+Nframes= numel(FRAMES);
+Nrows = ceil(Nframes/4);
+Ncols = ceil(Nframes/Nrows);
%
TNAME = [];
-FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 numel(FRAMES) 1])
+FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows])
for i_ = 1:numel(FRAMES)
- subplot(1,numel(FRAMES),i_); TNAME = [TNAME,'_',sprintf('%.0f',DATA.Ts3D(FRAMES(i_)))];
+ subplot(Nrows,Ncols,i_); TNAME = [TNAME,'_',sprintf('%.0f',DATA.Ts3D(FRAMES(i_)))];
if ~strcmp(OPTIONS.PLAN,'kxky')
scale = max(max(abs(toplot.FIELD(:,:,FRAMES(i_))))); % Scaling to normalize
else
scale = 1;
end
- pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale); set(pclr, 'edgecolor','none');pbaspect(toplot.ASPECT)
-% if ~strcmp(OPTIONS.PLAN,'kxky')
-% caxis([-1,1]);
-% colormap(bluewhitered);
-% end
+ pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale); set(pclr, 'edgecolor','none');
+ if OPTIONS.AXISEQUAL
+ pbaspect(toplot.ASPECT)
+ end
+ if ~strcmp(OPTIONS.PLAN,'kxky')
+ caxis([-1,1]);
+ colormap(bluewhitered);
+ end
xlabel(toplot.XNAME); ylabel(toplot.YNAME);
% if i_ > 1; set(gca,'ytick',[]); end;
title([sprintf('$t c_s/R=%.0f$',DATA.Ts3D(FRAMES(i_))),', max = ',sprintf('%.1e',scale)]);
diff --git a/matlab/plots/ZF_spacetime.m b/matlab/plots/ZF_spacetime.m
new file mode 100644
index 0000000000000000000000000000000000000000..b40deed90aa2ed66dbdf7aa5f0ff7761ccf7bec4
--- /dev/null
+++ b/matlab/plots/ZF_spacetime.m
@@ -0,0 +1,69 @@
+function [FIGURE] = ZF_spacetime(DATA, TAVG_0, TAVG_1)
+%Compute steady radial transport
+tend = TAVG_1; tstart = TAVG_0;
+[~,its3D] = min(abs(DATA.Ts3D-tstart));
+[~,ite3D] = min(abs(DATA.Ts3D-tend));
+nx = DATA.Nx; ny = DATA.Ny; nz = DATA.Nz; nkx = DATA.Nkx; nky = DATA.Nky;
+
+FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_spacetime','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1200, 600])
+
+%% radial shear radial profiles
+ % computation
+ Ns3D = numel(DATA.Ts3D);
+ [KY, KX] = meshgrid(DATA.ky, DATA.kx);
+ plt = @(x) mean(x(:,:,1,:),2);
+subplot(311)
+ phi = zeros(nx,ny,1,Ns3D);
+ for it = 1:numel(DATA.Ts3D)
+ % PERSONAL METHOD
+% phi(:,:,1,it) = real(fftshift(ifft2((DATA.PHI(:,:,1,it)),DATA.Nx,DATA.Ny)));
+ % GENE POST PROCESSING METHOD
+ phi(:,:,1,it) = ifourier_GENE(DATA.PHI(:,:,1,it),[DATA.Nx,DATA.Ny]);
+ end
+ clear spectrumKxKyZ spectrumKxKyZ
+ f2plot = phi; fname = '$\langle \phi\rangle_y$';
+ [TY,TX] = meshgrid(DATA.x,DATA.Ts3D);
+ pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
+ set(pclr, 'edgecolor','none');
+ legend(fname) %colorbar;
+ clim = max(max(abs(plt(f2plot(:,:,1,its3D:ite3D)))));
+ caxis(clim*[-1 1]);
+ cmap = bluewhitered(256);
+ colormap(cmap); colorbar;
+ ylabel('$x/\rho_s$');
+ title(DATA.param_title);
+
+subplot(312)
+ dxphi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+ for it = 1:numel(DATA.Ts3D)
+ dxphi(:,:,1,it) = ifourier_GENE(1i*KX.*DATA.PHI(:,:,1,it),[DATA.Nx,DATA.Ny]);
+ end
+ f2plot = dxphi; fname = '$\langle \partial_x\phi\rangle_y$';
+ [TY,TX] = meshgrid(DATA.x,DATA.Ts3D);
+ pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
+ set(pclr, 'edgecolor','none');
+ legend(fname) %colorbar;
+ clim = max(max(abs(plt(f2plot(:,:,1,its3D:ite3D)))));
+ caxis(clim*[-1 1]);
+ cmap = bluewhitered(256);
+ colormap(cmap); colorbar;
+ ylabel('$x/\rho_s$');
+
+subplot(313)
+ dx2phi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+ for it = 1:numel(DATA.Ts3D)
+ dx2phi(:,:,1,it) = ifourier_GENE(-KX.^2.*DATA.PHI(:,:,1,it),[DATA.Nx,DATA.Ny]);
+ end
+ f2plot = dx2phi;
+% fname = '$\Gamma_x(x)$';
+ fname = '$\langle \partial_x^2\phi\rangle_y$';
+ [TY,TX] = meshgrid(DATA.x,DATA.Ts3D);
+ pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
+ set(pclr, 'edgecolor','none');
+ legend(fname) %colorbar;
+ clim = max(max(abs(plt(f2plot(:,:,1,its3D:ite3D)))));
+ caxis(clim*[-1 1]);
+ cmap = bluewhitered(256);
+ colormap(cmap); colorbar;
+ xlabel('$t c_s/R$'), ylabel('$x/\rho_s$');
+end
\ No newline at end of file
diff --git a/matlab/plots/plot_radial_transport_and_shear.m b/matlab/plots/plot_radial_transport_and_shear.m
index 597952e0bfa91d795fff41c55300430aab81b9b7..97ca17ab42259b5d7f4087b74f426c9d5fa118b1 100644
--- a/matlab/plots/plot_radial_transport_and_shear.m
+++ b/matlab/plots/plot_radial_transport_and_shear.m
@@ -1,4 +1,4 @@
-function [FIGURE] = plot_radial_transport_and_shear(DATA, TAVG_0, TAVG_1)
+function [FIGURE] = plot_radial_transport_and_shear(DATA, TAVG_0, TAVG_1,stfname)
%Compute steady radial transport
tend = TAVG_1; tstart = TAVG_0;
[~,its0D] = min(abs(DATA.Ts0D-tstart));
@@ -18,26 +18,24 @@ FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; se
title(['$\nu_{',DATA.CONAME,'}=$', num2str(DATA.NU), ', $\kappa_N=$',num2str(DATA.K_N),...
', $L=',num2str(DATA.L),'$, $N=',num2str(DATA.Nx),'$, $(P,J)=(',num2str(DATA.PMAXI),',',num2str(DATA.JMAXI),')$,',...
' $\mu_{hd}=$',num2str(DATA.MU),', $\Gamma^{\infty} \approx $',num2str(gamma_infty_avg)]);
- %% zonal vs nonzonal energies for phi(t) and shear radial profile
-
+ %% zonal vs nonzonal energies for phi(t)
% computation
+ [~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(:,1,1,:))),2)));
Ns3D = numel(DATA.Ts3D);
[KY, KX] = meshgrid(DATA.ky, DATA.kx);
- Ephi_Z = zeros(1,Ns3D);
- Ephi_NZ_kgt0 = zeros(1,Ns3D);
- Ephi_NZ_kgt1 = zeros(1,Ns3D);
- Ephi_NZ_kgt2 = zeros(1,Ns3D);
- high_k_phi = zeros(1,Ns3D);
- dxphi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
-% dx2phi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+% Ephi_Z = zeros(1,Ns3D);
+% Ephi_NZ_kgt0 = zeros(1,Ns3D);
+% Ephi_NZ_kgt1 = zeros(1,Ns3D);
+% Ephi_NZ_kgt2 = zeros(1,Ns3D);
+ E_Zmode_SK = zeros(1,Ns3D);
+ E_NZmode_SK = zeros(1,Ns3D);
for it = 1:numel(DATA.Ts3D)
- [amp,ikzf] = max(abs((KX~=0).*DATA.PHI(:,1,1,it)));
- Ephi_NZ_kgt0(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>0.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2))));
- Ephi_NZ_kgt1(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>1.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2))));
- Ephi_NZ_kgt2(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>2.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2))));
- Ephi_Z(it) = squeeze(sum(sum(((KX~=0).*(KY==0).*(KX.^2).*abs(DATA.PHI(:,:,1,it)).^2))));
- dxphi(:,:,1,it) = real(fftshift(ifft2(-1i*KX.*(DATA.PHI(:,:,1,it)),DATA.Nx,DATA.Ny)));
-% dx2phi(:,:,1,it) = real(fftshift(ifft2(-KX.^2.*(DATA.PHI(:,:,1,it)),DATA.Nx,DATA.Ny)));
+% Ephi_NZ_kgt0(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>0.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2))));
+% Ephi_NZ_kgt1(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>1.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2))));
+% Ephi_NZ_kgt2(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>2.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2))));
+% Ephi_Z(it) = squeeze(sum(sum(((KX~=0).*(KY==0).*(KX.^2).*abs(DATA.PHI(:,:,1,it)).^2))));
+ E_Zmode_SK(it) = squeeze(DATA.ky(ikzf).^2.*abs(DATA.PHI(ikzf,1,1,it))^2);
+ E_NZmode_SK(it) = squeeze(sum(sum(((1+KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2.*(KY~=0)))));
end
[~,its3D] = min(abs(DATA.Ts3D-tstart));
[~,ite3D] = min(abs(DATA.Ts3D-tend));
@@ -45,21 +43,36 @@ FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; se
subplot(312)
it0 = 1; itend = Ns3D;
trange = it0:itend;
- pltx = @(x) x;%-x(1);
- plty = @(x) x./max((x(its3D:ite3D)));
+ plt1 = @(x) x;%-x(1);
+ plt2 = @(x) x./max((x(its3D:ite3D)));
% plty = @(x) x(500:end)./max(squeeze(x(500:end)));
yyaxis left
- plot(pltx(DATA.Ts3D(trange)),plty(Ephi_Z(trange)),'DisplayName',['Zonal, ',DATA.CONAME]); hold on;
+ plot(plt1(DATA.Ts3D(trange)),plt2(E_Zmode_SK(trange)),'DisplayName','$k_{zf}^2|\phi_{kzf}|^2$');
ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
yyaxis right
- plot(pltx(DATA.Ts3D(trange)),plty(Ephi_NZ_kgt0(trange)),'DisplayName','$k_p>0$');
+ plot(plt1(DATA.Ts3D(trange)),plt2(E_NZmode_SK(trange)),'DisplayName','$(1+k^2)|\phi^k|^2$');
xlim([DATA.Ts3D(it0), DATA.Ts3D(itend)]);
ylim([-0.1, 1.5]); ylabel('$E_{NZ}$')
xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
+ %% radial shear radial profile
+ % computation
+ Ns3D = numel(DATA.Ts3D);
+ [~, KX] = meshgrid(DATA.ky, DATA.kx);
+% dxphi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+% dx2phi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+ Gx = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+ for it = 1:numel(DATA.Ts3D)
+% dxphi(:,:,1,it) = real(fftshift(ifft2(-1i*KX.*(DATA.PHI(:,:,1,it)),DATA.Nx,DATA.Ny)));
+% dx2phi(:,:,1,it) = real(fftshift(ifft2(-KX.^2.*(DATA.PHI(:,:,1,it)),DATA.Nx,DATA.Ny)));
+ Gx(:,:,1,it) = real(fftshift(ifft2(-1i*KY.*(DATA.PHI(:,:,1,it)),DATA.Nx,DATA.Ny)))...
+ .*real(fftshift(ifft2(DATA.DENS_I(:,:,1,it),DATA.Nx,DATA.Ny)));
+ end
+ clim = max(max(abs(Gx)));
subplot(313)
[TY,TX] = meshgrid(DATA.x,DATA.Ts3D);
- pclr = pcolor(TX,TY,squeeze((mean(dxphi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_x\phi\rangle_y$') %colorbar;
+% pclr = pcolor(TX,TY,squeeze((mean(dxphi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_x\phi\rangle_y$') %colorbar;
% pclr = pcolor(TX,TY,squeeze((mean(dx2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_x^2\phi\rangle_y$') %colorbar;
- caxis([-3 3]);
+ pclr = pcolor(TX,TY,squeeze((mean(Gx(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_y \phi n_i\rangle_y$') %colorbar;
+ caxis(clim*[-1 1]);
xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); colormap(bluewhitered(256))
end
\ No newline at end of file
diff --git a/matlab/plots/plot_radial_transport_and_spacetime.m b/matlab/plots/plot_radial_transport_and_spacetime.m
new file mode 100644
index 0000000000000000000000000000000000000000..e6039a696796b110a20f78221903092d9c1c1288
--- /dev/null
+++ b/matlab/plots/plot_radial_transport_and_spacetime.m
@@ -0,0 +1,107 @@
+function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stfname)
+ %Compute steady radial transport
+ tend = TAVG_1; tstart = TAVG_0;
+ [~,its0D] = min(abs(DATA.Ts0D-tstart));
+ [~,ite0D] = min(abs(DATA.Ts0D-tend));
+ SCALE = (1/DATA.Nx/DATA.Ny)^2;
+ gamma_infty_avg = mean(DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
+ gamma_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
+ [~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(:,1,1,:))),2)));
+ Ns3D = numel(DATA.Ts3D);
+ [KY, KX] = meshgrid(DATA.ky, DATA.kx);
+
+ %% computations
+
+ % Compute Gamma from ifft matlab
+ Gx = zeros(DATA.Nx,DATA.Ny,numel(DATA.Ts3D));
+ for it = 1:numel(DATA.Ts3D)
+ for iz = 1:DATA.Nz
+ Gx(:,:,it) = Gx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)),[DATA.Nx,DATA.Ny])...
+ .*ifourier_GENE(DATA.DENS_I(:,:,iz,it),[DATA.Nx,DATA.Ny]);
+ end
+ Gx(:,:,it) = Gx(:,:,it)/DATA.Nz;
+ end
+ Gamma_t_mtlb = squeeze(mean(mean(Gx,1),2));
+ % zonal vs nonzonal energies for phi(t)
+
+ E_Zmode_SK = zeros(1,Ns3D);
+ E_NZmode_SK = zeros(1,Ns3D);
+ for it = 1:numel(DATA.Ts3D)
+ E_Zmode_SK(it) = squeeze(DATA.ky(ikzf).^2.*abs(DATA.PHI(ikzf,1,1,it))^2);
+ E_NZmode_SK(it) = squeeze(sum(sum(((1+KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2.*(KY~=0)))));
+ end
+ [~,its3D] = min(abs(DATA.Ts3D-tstart));
+ [~,ite3D] = min(abs(DATA.Ts3D-tend));
+
+%% Figure
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1200, 600])
+ subplot(311)
+% yyaxis left
+ plot(DATA.Ts0D,DATA.PGAMMA_RI*SCALE,'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
+ plot(DATA.Ts3D,Gamma_t_mtlb,'DisplayName','matlab comp.'); hold on;
+ plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',...
+ 'DisplayName',['$\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]);
+ grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_x$')
+ ylim([0,5*abs(gamma_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
+ title([DATA.param_title,', $\Gamma^{\infty} \approx $',num2str(gamma_infty_avg)]);
+ % plot
+ subplot(312)
+ it0 = 1; itend = Ns3D;
+ trange = it0:itend;
+ plt1 = @(x) x;%-x(1);
+ plt2 = @(x) x./max((x(its3D:ite3D)));
+% plty = @(x) x(500:end)./max(squeeze(x(500:end)));
+ yyaxis left
+ plot(plt1(DATA.Ts3D(trange)),plt2(E_Zmode_SK(trange)),'DisplayName','$k_{zf}^2|\phi_{kzf}|^2$');
+ ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
+ yyaxis right
+ plot(plt1(DATA.Ts3D(trange)),plt2(E_NZmode_SK(trange)),'DisplayName','$(1+k^2)|\phi^k|^2$');
+ xlim([DATA.Ts3D(it0), DATA.Ts3D(itend)]);
+ ylim([-0.1, 1.5]); ylabel('$E_{NZ}$')
+ xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
+ %% radial shear radial profile
+ % computation
+ Ns3D = numel(DATA.Ts3D);
+ [KY, KX] = meshgrid(DATA.ky, DATA.kx);
+ plt = @(x) mean(x(:,:,1,:),2);
+ switch stfname
+ case 'phi'
+ phi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+ for it = 1:numel(DATA.Ts3D)
+ phi(:,:,1,it) = ifourier_GENE(DATA.PHI(:,:,1,it),[DATA.Nx,DATA.Ny]);
+ end
+ f2plot = phi; fname = '$\phi$';
+ case 'v_y'
+ dxphi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+ for it = 1:numel(DATA.Ts3D)
+ dxphi(:,:,1,it) = ifourier_GENE(-1i*KX.*(DATA.PHI(:,:,1,it)),[DATA.Nx,DATA.Ny]);
+ end
+ f2plot = dxphi; fname = '$\langle \partial_x\phi\rangle_y$';
+ case 'v_x'
+ dyphi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+ for it = 1:numel(DATA.Ts3D)
+ dyphi(:,:,1,it) = ifourier_GENE(1i*KY.*(DATA.PHI(:,:,1,it)),[DATA.Nx,DATA.Ny]);
+ end
+ f2plot = dyphi; fname = '$\langle \partial_y\phi\rangle_y$';
+ case 'szf'
+ dx2phi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+ for it = 1:numel(DATA.Ts3D)
+ dx2phi(:,:,1,it) = ifourier_GENE(-KX.^2.*(DATA.PHI(:,:,1,it)),[DATA.Nx,DATA.Ny]);
+ end
+ f2plot = dx2phi; fname = '$\langle \partial_x^2\phi\rangle_y$';
+ end
+ clim = max(max(abs(plt(f2plot(:,:,1,its3D:ite3D)))));
+ subplot(313)
+ [TY,TX] = meshgrid(DATA.x,DATA.Ts3D);
+ pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
+ set(pclr, 'edgecolor','none');
+ legend(fname) %colorbar;
+ caxis(clim*[-1 1]);
+ cmap = bluewhitered(256);
+ colormap(cmap)
+ xlabel('$t c_s/R$'), ylabel('$x/\rho_s$');
+ if strcmp(stfname,'Gx')
+ subplot(311)
+ plot(DATA.Ts3D,squeeze(mean(plt(f2plot),1)));
+ end
+end
\ No newline at end of file
diff --git a/matlab/process_field.m b/matlab/process_field.m
index b2df3a5f5477454497307b3dd3270c62e46af163..376186a7c062ecea479485046544d2acdb1fd679 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -147,6 +147,69 @@ switch OPTIONS.NAME
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
+ case 'n_e'
+ NAME = 'ne';
+ if COMPDIM == 3
+ for it = FRAMES
+ tmp = squeeze(compr(DATA.DENS_E(:,:,:,it)));
+ FIELD(:,:,it) = squeeze(process(tmp));
+ end
+ else
+ if REALP
+ tmp = zeros(Nx,Ny,Nz);
+ else
+ tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ end
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,it)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
+ case 'k^2n_e'
+ NAME = 'k2ne';
+ [KY, KX] = meshgrid(DATA.ky, DATA.kx);
+ if COMPDIM == 3
+ for it = FRAMES
+ for iz = 1:DATA.Nz
+ tmp = squeeze(compr(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,it)));
+ end
+ FIELD(:,:,it) = squeeze(process(tmp));
+ end
+ else
+ if REALP
+ tmp = zeros(Nx,Ny,Nz);
+ else
+ tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ end
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = squeeze(process(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,it)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
+ case 'n_e^{NZ}'
+ NAME = 'neNZ';
+ if COMPDIM == 3
+ for it = FRAMES
+ tmp = squeeze(compr(DATA.DENS_E(:,:,:,it).*(KY~=0)));
+ FIELD(:,:,it) = squeeze(process(tmp));
+ end
+ else
+ if REALP
+ tmp = zeros(Nx,Ny,Nz);
+ else
+ tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ end
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,it).*(KY~=0)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
case 'n_i'
NAME = 'ni';
if COMPDIM == 3
@@ -207,58 +270,78 @@ switch OPTIONS.NAME
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
- case '\Gamma_y'
- NAME = 'Gy';
- [KY, ~] = meshgrid(DATA.ky, DATA.kx);
- Gx = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
+ case 'v_y'
+ NAME = 'vy';
+ [~, KX] = meshgrid(DATA.ky, DATA.kx);
+ vE = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
for it = FRAMES % Compute the 3D real transport for each frame
for iz = 1:DATA.Nz
- Gx(:,:,iz,it) = real(ifft2(-1i*KY.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny))...
- .*real(ifft2(DATA.DENS_I(:,:,iz,it),DATA.Nx,DATA.Ny));
+ vE(:,:,iz,it) = real(ifft2(-1i*KX.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny));
end
end
if REALP % plot in real space
for it = FRAMES
- FIELD(:,:,it) = squeeze(compr(Gx(:,:,:,it)));
+ FIELD(:,:,it) = squeeze(compr(vE(:,:,:,it)));
end
else % Plot spectrum
process = @(x) abs(fftshift(x,2));
- shift_x = @(x) fftshift(x,2);
- shift_y = @(x) fftshift(x,2);
tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
for it = FRAMES
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = process(squeeze(fft2(Gx(:,:,iz,it),DATA.Nkx,DATA.Nky)));
+ tmp(:,:,iz) = squeeze(process(-1i*KX.*(DATA.PHI(:,:,iz,it))));
end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case '\Gamma_x'
- NAME = 'Gx';
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
+ case 'v_x'
+ NAME = 'vx';
[KY, ~] = meshgrid(DATA.ky, DATA.kx);
- Gx = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
+ vE = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
for it = FRAMES % Compute the 3D real transport for each frame
for iz = 1:DATA.Nz
- Gx(:,:,iz,it) = real((ifft2(-1i*KY.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny)))...
- .*real((ifft2(DATA.DENS_I(:,:,iz,it),DATA.Nx,DATA.Ny)));
+ vE(:,:,iz,it) = real(ifft2(-1i*KY.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny));
end
end
if REALP % plot in real space
for it = FRAMES
- FIELD(:,:,it) = squeeze(compr(Gx(:,:,:,it)));
+ FIELD(:,:,it) = squeeze(compr(vE(:,:,:,it)));
end
else % Plot spectrum
process = @(x) abs(fftshift(x,2));
- shift_x = @(x) fftshift(x,2);
- shift_y = @(x) fftshift(x,2);
- tmp = zeros(DATA.Nx,DATA.Ny,Nz);
+ tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
for it = FRAMES
for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = process(squeeze(fft2(Gx(:,:,iz,it),DATA.Nx,DATA.Ny)));
+ tmp(:,:,iz) = squeeze(process(-1i*KY.*(DATA.PHI(:,:,iz,it))));
end
- FIELD(:,:,it) = squeeze(compr(tmp(1:DATA.Nkx,1:DATA.Nky,:)));
- end
- end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
+ case '\Gamma_x'
+ NAME = 'Gx';
+ [KY, ~] = meshgrid(DATA.ky, DATA.kx);
+ Gx = zeros(DATA.Nx,DATA.Ny,DATA.Nz,numel(FRAMES));
+ for it = FRAMES % Compute the 3D real transport for each frame
+ for iz = 1:DATA.Nz
+ Gx(:,:,iz,it) = real((ifft2(-1i*KY.*(DATA.PHI(:,:,iz,it)),DATA.Nx,DATA.Ny)))...
+ .*real((ifft2(DATA.DENS_I(:,:,iz,it),DATA.Nx,DATA.Ny)));
+ end
+ end
+ if REALP % plot in real space
+ for it = FRAMES
+ FIELD(:,:,it) = squeeze(compr(Gx(:,:,:,it)));
+ end
+ else % Plot spectrum
+ process = @(x) abs(fftshift(x,2));
+ shift_x = @(x) fftshift(x,2);
+ shift_y = @(x) fftshift(x,2);
+ tmp = zeros(DATA.Nx,DATA.Ny,Nz);
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = process(squeeze(fft2(Gx(:,:,iz,it),DATA.Nx,DATA.Ny)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp(1:DATA.Nkx,1:DATA.Nky,:)));
+ end
+ end
end
TOPLOT.FIELD = FIELD;
TOPLOT.X = shift_x(X);
diff --git a/matlab/setup.m b/matlab/setup.m
index f696387e5cf5109889e5d32c83f20bb6fc288d57..5ec853a92ab44795775c5b6e484dc81f10cf29cb 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -17,10 +17,12 @@ if SG; GRID.SG = '.true.'; else; GRID.SG = '.false.';end;
% Model parameters
MODEL.CLOS = CLOS;
MODEL.NL_CLOS = NL_CLOS;
-MODEL.LINEARITY = LINEARITY;
+MODEL.LINEARITY = ['''',LINEARITY,''''];
MODEL.KIN_E = KIN_E;
if KIN_E; MODEL.KIN_E = '.true.'; else; MODEL.KIN_E = '.false.';end;
-MODEL.mu = MU;
+MODEL.mu_x = MU_X;
+MODEL.mu_y = MU_Y;
+MODEL.mu_z = 0;
MODEL.mu_p = MU_P;
MODEL.mu_j = MU_J;
MODEL.nu = NU; % hyper diffusive coefficient nu for HW
@@ -42,7 +44,7 @@ MODEL.GradB = GRADB; % Magnetic gradient
MODEL.CurvB = CURVB; % Magnetic curvature
MODEL.lambdaD = LAMBDAD;
% Collision parameters
-COLL.collision_model = CO;
+COLL.collision_model = ['''',CO,''''];
if (GKCO); COLL.gyrokin_CO = '.true.'; else; COLL.gyrokin_CO = '.false.';end;
if (ABCO); COLL.interspecies = '.true.'; else; COLL.interspecies = '.false.';end;
COLL.mat_file = '''null''';
@@ -56,14 +58,11 @@ switch CO
end
% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4''';
-if (INIT_PHI); INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end;
-INITIAL.INIT_ZF = INIT_ZF;
-INITIAL.wipe_turb = WIPE_TURB;
-INITIAL.ACT_ON_MODES = ACT_ON_MODES;
-if (INIT_BLOB); INITIAL.init_blob = '.true.'; else; INITIAL.init_blob = '.false.';end;
-INITIAL.init_background = (INIT_ZF>0)*ZF_AMP + BCKGD0;
-INITIAL.init_noiselvl = NOISE0;
-INITIAL.iseed = 42;
+INITIAL.INIT_OPT = ['''',INIT_OPT,''''];
+INITIAL.ACT_ON_MODES = ['''',ACT_ON_MODES,''''];
+INITIAL.init_background = BCKGD0;
+INITIAL.init_noiselvl = NOISE0;
+INITIAL.iseed = 42;
% Naming and creating input file
CONAME = CO;
diff --git a/matlab/show_geometry.m b/matlab/show_geometry.m
index e57746d858498f5a68aaa27708195c04b8b7c78f..3794d8c8ec17a61b29391e0d000142b82381d9c9 100644
--- a/matlab/show_geometry.m
+++ b/matlab/show_geometry.m
@@ -12,22 +12,20 @@ if DATA.Nz > 1 % Torus flux tube geometry
x_tor = (R + a.*cos(p)) .* cos(t);
y_tor = (R + a.*cos(p)) .* sin(t);
z_tor = a.*sin(p);
+ DIMENSIONS = [50 50 1200 600];
else % Zpinch geometry
Nturns = 0.1; Tgeom = 0;
- a = 1; % Torus minor radius
+ a = 0.7; % Torus minor radius
R = 0; % Torus major radius
q = 0;
theta = linspace(-Nturns*pi, Nturns*pi, 100); % Poloidal angle
- phi = linspace(-0.1, 0.1, 3); % cylinder height
+ phi = linspace(-0.1, 0.1, 5); % cylinder height
[t, p] = meshgrid(phi, theta);
x_tor = a.*cos(p);
z_tor = a.*sin(p);
y_tor = t;
+ DIMENSIONS = [50 50 numel(OPTIONS.TIME)*400 500];
end
-FIGURE.fig = figure; set(gcf, 'Position', [50 50 1200 600])
-torus=surf(x_tor, y_tor, z_tor); hold on;alpha 1.0;%light('Position',[-1 1 1],'Style','local')
-set(torus,'edgecolor',[1 1 1]*0.8,'facecolor','none')
-xlabel('X');ylabel('Y');zlabel('Z');
% field line coordinates
Xfl = @(z) (R+a*cos(z)).*cos(q*z);
Yfl = @(z) (R+a*cos(z)).*sin(q*z);
@@ -49,16 +47,10 @@ yY = @(z) bZ(z).*xX(z)-bX(z).*xZ(z)./sqrt((bY(z).*xZ(z)-bZ(z).*xY(z)).^2+(bZ(z)
yZ = @(z) bX(z).*xY(z)-bY(z).*xX(z)./sqrt((bY(z).*xZ(z)-bZ(z).*xY(z)).^2+(bZ(z).*xX(z)-bX(z).*xZ(z)).^2+(bX(z).*xY(z)-bY(z).*xX(z)).^2);
yvec= @(z) [yX(z); yY(z); yZ(z)];
% Plot high res field line
-theta = linspace(-Nturns*pi, Nturns*pi, 512) ; % Poloidal angle
-plot3(Xfl(theta),Yfl(theta),Zfl(theta))
% Planes plot
theta = DATA.z ; % Poloidal angle
-plot3(Xfl(theta),Yfl(theta),Zfl(theta),'ok')
% Plot x,y,bvec at these points
-scale = 0.15;
-quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*xX(theta),scale*xY(theta),scale*xZ(theta),0,'r');
-quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*yX(theta),scale*yY(theta),scale*yZ(theta),0,'g');
-quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*bX(theta),scale*bY(theta),scale*bZ(theta),0,'b');
+scale = 0.10;
%% Plot plane result
OPTIONS.POLARPLOT = 0;
@@ -66,31 +58,52 @@ OPTIONS.PLAN = 'xy';
r_o_R = DATA.rho_o_R;
[Y,X] = meshgrid(r_o_R*DATA.y,r_o_R*DATA.x);
max_ = 0;
-FIELDS = zeros(DATA.Nx,DATA.Ny,numel(OPTIONS.PLANES));
+FIELDS = zeros(DATA.Nx,DATA.Ny,DATA.Nz);
+FIGURE.fig = figure; FIGURE.FIGNAME = ['geometry','_',DATA.PARAMS]; set(gcf, 'Position', DIMENSIONS)
for it_ = 1:numel(OPTIONS.TIME)
-[~,it] = min(abs(OPTIONS.TIME(it_)-DATA.Ts3D));
-for iz = OPTIONS.PLANES
- OPTIONS.COMP = iz;
- toplot = process_field(DATA,OPTIONS);
- FIELDS(:,:,iz) = toplot.FIELD(:,:,it);
- tmp = max(max(abs(FIELDS(:,:,iz))));
- if (tmp > max_) max_ = tmp;
-end
-for iz = OPTIONS.PLANES
- z_ = DATA.z(iz);
- Xp = Xfl(z_) + X*xX(z_) + Y*yX(z_);
- Yp = Yfl(z_) + X*xY(z_) + Y*yY(z_);
- Zp = Zfl(z_) + X*xZ(z_) + Y*yZ(z_);
- s=surface(Xp,Yp,Zp,FIELDS(:,:,iz)/max_);
- s.EdgeColor = 'none';
- colormap(bluewhitered);
- caxis([-1,1]);
-end
-end
-%%
-axis equal
-view([1,-2,1])
+subplot(1,numel(OPTIONS.TIME),it_)
+ %plot magnetic geometry
+ 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.8,'facecolor','none')
+ %plot field line
+ theta = linspace(-Nturns*pi, Nturns*pi, 512) ; % Poloidal angle
+ plot3(Xfl(theta),Yfl(theta),Zfl(theta))
+ %plot vector basis
+ theta = DATA.z ; % Poloidal angle
+ plot3(Xfl(theta),Yfl(theta),Zfl(theta),'ok'); hold on;
+ quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*xX(theta),scale*xY(theta),scale*xZ(theta),0,'r');
+ quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*yX(theta),scale*yY(theta),scale*yZ(theta),0,'g');
+ quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*bX(theta),scale*bY(theta),scale*bZ(theta),0,'b');
+ xlabel('X');ylabel('Y');zlabel('Z');
+ %Plot time dependent data
+ [~,it] = min(abs(OPTIONS.TIME(it_)-DATA.Ts3D));
+ for iz = OPTIONS.PLANES
+ OPTIONS.COMP = iz;
+ toplot = process_field(DATA,OPTIONS);
+ FIELDS(:,:,iz) = toplot.FIELD(:,:,it);
+ tmp = max(max(abs(FIELDS(:,:,iz))));
+ if (tmp > max_) max_ = tmp;
+ end
+ for iz = OPTIONS.PLANES
+ z_ = DATA.z(iz);
+ Xp = Xfl(z_) + X*xX(z_) + Y*yX(z_);
+ Yp = Yfl(z_) + X*xY(z_) + Y*yY(z_);
+ Zp = Zfl(z_) + X*xZ(z_) + Y*yZ(z_);
+ s=surface(Xp,Yp,Zp,FIELDS(:,:,iz)/max_);
+ s.EdgeColor = 'none';
+ colormap(bluewhitered);
+% caxis([-1,1]);
+ end
+ if DATA.Nz == 1
+ xlim([0.65 0.75]);
+ ylim([-0.1 0.1]);
+ zlim([-0.2 0.2]);
+ end
+ end
+ %%
+ % axis equal
+ view([1,-2,1])
end
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index e7901dd410458c0c333331425e12e0bba7b19e8a..574b6d02035fc7e72a789a5c131e2ef24bb2225b 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -50,7 +50,9 @@ fprintf(fid,[' CLOS = ', num2str(MODEL.CLOS),'\n']);
fprintf(fid,[' NL_CLOS = ', num2str(MODEL.NL_CLOS),'\n']);
fprintf(fid,[' LINEARITY = ', MODEL.LINEARITY,'\n']);
fprintf(fid,[' KIN_E = ', MODEL.KIN_E,'\n']);
-fprintf(fid,[' mu = ', num2str(MODEL.mu),'\n']);
+fprintf(fid,[' mu_x = ', num2str(MODEL.mu_x),'\n']);
+fprintf(fid,[' mu_y = ', num2str(MODEL.mu_y),'\n']);
+fprintf(fid,[' mu_z = ', num2str(MODEL.mu_z),'\n']);
fprintf(fid,[' mu_p = ', num2str(MODEL.mu_p),'\n']);
fprintf(fid,[' mu_j = ', num2str(MODEL.mu_j),'\n']);
fprintf(fid,[' nu = ', num2str(MODEL.nu),'\n']);
@@ -77,11 +79,8 @@ fprintf(fid,'/\n');
fprintf(fid,'&INITIAL_CON\n');
-fprintf(fid,[' INIT_NOISY_PHI = ', INITIAL.init_noisy_phi,'\n']);
-fprintf(fid,[' INIT_ZF = ', num2str(INITIAL.INIT_ZF),'\n']);
-fprintf(fid,[' ACT_ON_MODES = ', INITIAL.ACT_ON_MODES,'\n']);
-fprintf(fid,[' WIPE_TURB = ', num2str(INITIAL.wipe_turb),'\n']);
-fprintf(fid,[' INIT_BLOB = ', INITIAL.init_blob,'\n']);
+fprintf(fid,[' INIT_OPT = ', INITIAL.INIT_OPT,'\n']);
+fprintf(fid,[' ACT_ON_MODES = ', INITIAL.ACT_ON_MODES,'\n']);
fprintf(fid,[' init_background = ', num2str(INITIAL.init_background),'\n']);
fprintf(fid,[' init_noiselvl = ', num2str(INITIAL.init_noiselvl),'\n']);
fprintf(fid,[' iseed = ', num2str(INITIAL.iseed),'\n']);
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index 835c9221c0e2d611b5c618c926b180968afed989..9bdcd861187e940466122d53b61650a74df73bfb 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -49,6 +49,9 @@ CONTAINS
interspecies = .false.
CASE ('LD') ! Landau
cosolver_coll = .true.
+ CASE DEFAULT
+ cosolver_coll = .false.
+ interspecies = .false.
END SELECT
END SUBROUTINE collision_readinputs
@@ -72,18 +75,20 @@ CONTAINS
SUBROUTINE compute_TColl
USE basic
+ USE model, ONLY : nu
IMPLICIT NONE
! Execution time start
CALL cpu_time(t0_coll)
-
- SELECT CASE(collision_model)
- CASE ('LB')
- CALL compute_lenard_bernstein
- CASE ('DG')
- CALL compute_dougherty
- CASE DEFAULT
- CALL compute_cosolver_coll
- END SELECT
+ IF (nu .NE. 0) THEN
+ SELECT CASE(collision_model)
+ CASE ('LB')
+ CALL compute_lenard_bernstein
+ CASE ('DG')
+ CALL compute_dougherty
+ CASE DEFAULT
+ CALL compute_cosolver_coll
+ END SELECT
+ ENDIF
! Execution time end
CALL cpu_time(t1_coll)
diff --git a/src/inital.F90 b/src/inital.F90
index 880143c8187543a9737995af62c3cb28a31f994e..06a336b0f6ebecc7e11c6727d27a837ea70e7d6a 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -13,7 +13,7 @@ SUBROUTINE inital
USE closure
USE ghosts
USE restarts
- USE numerics, ONLY: wipe_turbulence, play_with_modes, save_EM_ZF_modes
+ USE numerics, ONLY: play_with_modes, save_EM_ZF_modes
USE processing, ONLY: compute_nadiab_moments
IMPLICIT NONE
@@ -27,17 +27,26 @@ SUBROUTINE inital
CALL poisson ! compute phi_0=phi(N_0)
! through initialization
ELSE
+ SELECT CASE (INIT_OPT)
! set phi with noise and set moments to 0
- IF (INIT_NOISY_PHI) THEN
+ CASE ('phi')
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
CALL init_phi
! set moments_00 (GC density) with noise and compute phi afterwards
- ELSE
+ CASE('mom00')
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy gyrocenter density'
CALL init_gyrodens ! init only gyrocenter density
! CALL init_moments ! init all moments randomly (unadvised)
CALL poisson ! get phi_0 = phi(N_0)
- ENDIF
+ CASE('allmom')
+ IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy moments'
+ CALL init_moments ! init all moments
+ CALL poisson ! get phi_0 = phi(N_0)
+ CASE('blob')
+ IF (my_id .EQ. 0) WRITE(*,*) '--init a blob'
+ CALL initialize_blob
+ CALL poisson ! get phi_0 = phi(N_0)
+ END SELECT
ENDIF
! Save zonal and entropy modes
@@ -46,12 +55,6 @@ SUBROUTINE inital
! Freeze/Wipe some selected modes (entropy,zonal,turbulent)
CALL play_with_modes
- ! Option for initializing a gaussian blob on the zonal profile
- IF ( INIT_BLOB ) THEN
- IF (my_id .EQ. 0) WRITE(*,*) '--init a blob'
- CALL initialize_blob
- ENDIF
-
IF (my_id .EQ. 0) WRITE(*,*) 'Apply closure'
CALL apply_closure_model
@@ -338,6 +341,16 @@ SUBROUTINE initialize_blob
ENDIF
ENDDO
ENDDO
+ DO ip=ips_e,ipe_e
+ DO ij=ijs_e,ije_e
+ IF( (iky .NE. iky_0) .AND. (ip .EQ. 1) .AND. (ij .EQ. 1)) THEN
+ moments_e( ip,ij,ikx,iky,iz, :) = moments_e( ip,ij,ikx,iky,iz, :) &
+ + gain*sigma/SQRT2 * exp(-(kx**2+ky**2)*sigma**2/4._dp) &
+ * AA_x(ikx)*AA_y(iky)!&
+ ! * exp(sigmai2_taui_o2*(kx**2+ky**2))
+ ENDIF
+ ENDDO
+ ENDDO
ENDDO
ENDDO
ENDDO
diff --git a/src/initial_par_mod.F90 b/src/initial_par_mod.F90
index fc72d99e9a7afea02f9955391653819a100bfc8c..a9d9714f62728f5902bfc7f88c3e4c2f0c3d6385 100644
--- a/src/initial_par_mod.F90
+++ b/src/initial_par_mod.F90
@@ -5,17 +5,13 @@ MODULE initial_par
IMPLICIT NONE
PRIVATE
- ! Initialization through a noisy phi
- LOGICAL, PUBLIC, PROTECTED :: INIT_NOISY_PHI = .false.
+ ! Initialization option (phi/mom00/allmom/blob)
+ CHARACTER(len=32), PUBLIC, PROTECTED :: INIT_OPT = 'phi'
! Initialization through a zonal flow phi
INTEGER, PUBLIC, PROTECTED :: INIT_ZF = 0
REAL(DP), PUBLIC, PROTECTED :: ZF_AMP = 1E+3_dp
! Act on modes artificially (keep/wipe, zonal, non zonal, entropy mode etc.)
CHARACTER(len=32), PUBLIC, PROTECTED :: ACT_ON_MODES = 'nothing'
- ! Wipe turbulence in the restart (=1) or at each step (=2)
- INTEGER, PUBLIC, PROTECTED :: WIPE_TURB = 0
- ! Init a Gaussian blob density in the middle
- LOGICAL, PUBLIC, PROTECTED :: INIT_BLOB = .false.
! Initial background level
REAL(dp), PUBLIC, PROTECTED :: init_background=0._dp
! Initial noise amplitude
@@ -35,11 +31,8 @@ CONTAINS
USE prec_const
IMPLICIT NONE
- NAMELIST /INITIAL_CON/ INIT_NOISY_PHI
- NAMELIST /INITIAL_CON/ INIT_ZF
+ NAMELIST /INITIAL_CON/ INIT_OPT
NAMELIST /INITIAL_CON/ ACT_ON_MODES
- NAMELIST /INITIAL_CON/ WIPE_TURB
- NAMELIST /INITIAL_CON/ INIT_BLOB
NAMELIST /INITIAL_CON/ init_background
NAMELIST /INITIAL_CON/ init_noiselvl
NAMELIST /INITIAL_CON/ iseed
@@ -59,9 +52,7 @@ CONTAINS
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
- CALL attach(fidres, TRIM(str), "INIT_NOISY_PHI", INIT_NOISY_PHI)
-
- CALL attach(fidres, TRIM(str), "INIT_ZF", INIT_ZF)
+ CALL attach(fidres, TRIM(str), "INIT_OPT", INIT_OPT)
CALL attach(fidres, TRIM(str), "init_background", init_background)
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 2f56451f96761098c4a91da9c35d8558e963c160..2e1f4d9d055da973e9dd5476fb7bb10b3f5b4e8a 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -11,7 +11,9 @@ MODULE model
CHARACTER(len=16), &
PUBLIC, PROTECTED ::LINEARITY= 'linear' ! To turn on non linear bracket term
LOGICAL, PUBLIC, PROTECTED :: KIN_E = .true. ! Simulate kinetic electron (adiabatic otherwise)
- REAL(dp), PUBLIC, PROTECTED :: mu = 0._dp ! spatial Hyperdiffusivity coefficient (for num. stability)
+ REAL(dp), PUBLIC, PROTECTED :: mu_x = 0._dp ! spatial x-Hyperdiffusivity coefficient (for num. stability)
+ REAL(dp), PUBLIC, PROTECTED :: mu_y = 0._dp ! spatial y-Hyperdiffusivity coefficient (for num. stability)
+ REAL(dp), PUBLIC, PROTECTED :: mu_z = 0._dp ! spatial z-Hyperdiffusivity coefficient (for num. stability)
REAL(dp), PUBLIC, PROTECTED :: mu_p = 0._dp ! kinetic para hyperdiffusivity coefficient (for num. stability)
REAL(dp), PUBLIC, PROTECTED :: mu_j = 0._dp ! kinetic perp hyperdiffusivity coefficient (for num. stability)
REAL(dp), PUBLIC, PROTECTED :: nu = 1._dp ! Collision frequency
@@ -54,8 +56,10 @@ CONTAINS
USE prec_const
IMPLICIT NONE
- NAMELIST /MODEL_PAR/ CLOS, NL_CLOS, KERN, LINEARITY, KIN_E, mu, mu_p, mu_j, nu, tau_e, tau_i, sigma_e, sigma_i, &
- q_e, q_i, K_n, K_T, K_E, GradB, CurvB, lambdaD
+ NAMELIST /MODEL_PAR/ CLOS, NL_CLOS, KERN, LINEARITY, KIN_E, &
+ mu_x, mu_y, mu_z, mu_p, mu_j, nu,&
+ tau_e, tau_i, sigma_e, sigma_i, q_e, q_i,&
+ K_n, K_T, K_E, GradB, CurvB, lambdaD
READ(lu_in,model_par)
@@ -104,7 +108,10 @@ CONTAINS
CALL attach(fidres, TRIM(str), "LINEARITY", LINEARITY)
CALL attach(fidres, TRIM(str), "KIN_E", KIN_E)
CALL attach(fidres, TRIM(str), "nu", nu)
- CALL attach(fidres, TRIM(str), "mu", mu)
+ CALL attach(fidres, TRIM(str), "mu", 0)
+ CALL attach(fidres, TRIM(str), "mu_x", mu_x)
+ CALL attach(fidres, TRIM(str), "mu_y", mu_y)
+ CALL attach(fidres, TRIM(str), "mu_z", mu_z)
CALL attach(fidres, TRIM(str), "tau_e", tau_e)
CALL attach(fidres, TRIM(str), "tau_i", tau_i)
CALL attach(fidres, TRIM(str), "sigma_e", sigma_e)
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 29f34e8035c0403523c186be8d0f737964bc070f..4e5ce055ac020b85a6451039ec47615c1d8d6118 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -121,7 +121,7 @@ SUBROUTINE moments_eq_rhs_e
! Electrostatic background gradients
- i_ky * K_E * moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
! Numerical hyperdiffusion (totally artificial, for stability purpose)
- - mu*kperp2**2 * moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
+ - (mu_x*kx**2 + mu_y*ky**2)*moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
! Collision term
+ TColl_e(ip,ij,ikx,iky,iz)
@@ -258,7 +258,7 @@ SUBROUTINE moments_eq_rhs_i
! Electrostatic background gradients
- i_ky * K_E * moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
! Numerical hyperdiffusion (totally artificial, for stability purpose)
- - mu*kperp2**2 * moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
+ - (mu_x*kx**2 + mu_y*ky**2)*moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
! Collision term
+ TColl_i(ip,ij,ikx,iky,iz)
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 1b3926173c113de07eaf802f64ae1eb35d8c3dbd..84c1b30a19409a5a14dba0f90b98d88dd1cd3395 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -7,7 +7,7 @@ MODULE numerics
implicit none
PUBLIC :: build_dnjs_table, evaluate_kernels, evaluate_poisson_op, compute_lin_coeff
- PUBLIC :: wipe_turbulence, play_with_modes, save_EM_ZF_modes
+ PUBLIC :: play_with_modes, save_EM_ZF_modes
CONTAINS
@@ -244,43 +244,6 @@ SUBROUTINE compute_lin_coeff
END SUBROUTINE compute_lin_coeff
-!******************************************************************************!
-!!!!!!! Remove all ky!=0 modes to conserve only zonal modes in a restart
-!******************************************************************************!
-SUBROUTINE wipe_turbulence
- USE fields
- USE grid
- IMPLICIT NONE
- DO ikx=ikxs,ikxe
- DO iky=ikys,ikye
- DO iz=izs,ize
- DO ip=ips_e,ipe_e
- DO ij=ijs_e,ije_e
- IF( iky .NE. iky_0) THEN
- moments_e( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_e( ip,ij,ikx,iky,iz, :)
- ELSE
- moments_e( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_e( ip,ij,ikx,iky,iz, :)
- ENDIF
- ENDDO
- ENDDO
- DO ip=ips_i,ipe_i
- DO ij=ijs_i,ije_i
- IF( iky .NE. iky_0) THEN
- moments_i( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_i( ip,ij,ikx,iky,iz, :)
- ELSE
- moments_i( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_i( ip,ij,ikx,iky,iz, :)
- ENDIF
- ENDDO
- ENDDO
- IF( iky .NE. iky_0) THEN
- phi(ikx,iky,iz) = 0e-3_dp*phi(ikx,iky,iz)
- ELSE
- phi(ikx,iky,iz) = 1e+0_dp*phi(ikx,iky,iz)
- ENDIF
- ENDDO
- ENDDO
- ENDDO
-END SUBROUTINE
!******************************************************************************!
!!!!!!! Routine that can artificially increase or wipe modes
!******************************************************************************!
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index 0bb98587c82da2e93cb14d1c7e1ecdd797d277f5..797267b88d2c35d6f70f3d551b4eb40fc1125dfb 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -1,6 +1,7 @@
MODULE restarts
USE basic
-USE futils, ONLY: openf, closef, getarr, getatt, isgroup, isdataset,getarrnd
+USE futils, ONLY: openf, closef, getarr, getatt, isgroup,&
+ isdataset, getarrnd, putarrnd
USE grid
USE fields
USE diagnostics_par
@@ -15,7 +16,7 @@ INTEGER :: pmaxe_cp, jmaxe_cp, pmaxi_cp, jmaxi_cp, n0
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_e_cp
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_i_cp
-PUBLIC :: load_moments
+PUBLIC :: load_moments!, write_restart
CONTAINS
diff --git a/src/srcinfo.h b/src/srcinfo.h
index 7cc5e9ad5cd3212593dd6b361c4cf1a77640ab8f..0d57eb336ae8322c3f325f0a5f265a715f5bee74 100644
--- a/src/srcinfo.h
+++ b/src/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='2df4eb7-dirty')
+parameter (VERSION='6d12f9c-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Tue Nov 23 10:32:57 CET 2021')
+parameter (EXECDATE='Tue Dec 14 10:47:37 CET 2021')
parameter (HOST ='spcpc606')
diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h
index 7cc5e9ad5cd3212593dd6b361c4cf1a77640ab8f..0d57eb336ae8322c3f325f0a5f265a715f5bee74 100644
--- a/src/srcinfo/srcinfo.h
+++ b/src/srcinfo/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='2df4eb7-dirty')
+parameter (VERSION='6d12f9c-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Tue Nov 23 10:32:57 CET 2021')
+parameter (EXECDATE='Tue Dec 14 10:47:37 CET 2021')
parameter (HOST ='spcpc606')
diff --git a/src/stepon.F90 b/src/stepon.F90
index a57b3e4bde6d2e6391f583d9b596dd80096d7d1a..f56880aee24408625d08627b5cb56cd9a6081752 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -6,13 +6,13 @@ SUBROUTINE stepon
USE closure
USE collision, ONLY : compute_TColl
USE fields, ONLY: moments_e, moments_i, phi
- USE initial_par, ONLY: ACT_ON_MODES, WIPE_TURB
+ USE initial_par, ONLY: ACT_ON_MODES
USE ghosts
USE grid
USE model, ONLY : LINEARITY, KIN_E
use prec_const
USE time_integration
- USE numerics, ONLY: play_with_modes, wipe_turbulence
+ USE numerics, ONLY: play_with_modes
USE processing, ONLY: compute_nadiab_moments
USE utility, ONLY: checkfield
@@ -47,8 +47,6 @@ SUBROUTINE stepon
CALL compute_Sapj
! Store or cancel/maintain zonal modes artificially
CALL play_with_modes
- ! Cancel non zonal modes artificially
- IF ( WIPE_TURB .EQ. 2) CALL wipe_turbulence
!- Check before next step
CALL checkfield_all()
IF( nlend ) EXIT ! exit do loop
diff --git a/wk/Dimits_shift_results.m b/wk/Dimits_shift_results.m
new file mode 100644
index 0000000000000000000000000000000000000000..5bbd02a31c811713a48a5e41f9076be9b0de41a3
--- /dev/null
+++ b/wk/Dimits_shift_results.m
@@ -0,0 +1,88 @@
+%% Nu = 1e-2
+KN = [ 1.5 1.6 1.7 1.8 1.9];
+Ginf_LD = [5.5e-2 1.3e-1 2.7e-1 4.6e-1 1.0e+0];
+conv_LD = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+
+Ginf_SG = [1.2e-3 2.8e-3 1.7e-2 8.9e-2 3.0e+0];
+conv_SG = [9.0e-4 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+
+
+Ginf_DG = [1.1e-4 5.6e-4 9.0e-4 6.5e-3 7.0e+0];
+conv_DG = [1.6e-4 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+
+figure
+semilogy(KN,Ginf_LD,'d--g','DisplayName','Landau'); hold on;
+
+
+semilogy(KN,Ginf_SG,'o-b','DisplayName','Sugama'); hold on;
+semilogy(KN,conv_SG,'x-k','DisplayName','conv'); hold on;
+
+
+semilogy(KN,Ginf_DG,'^-r','DisplayName','Dougherty'); hold on;
+semilogy(KN,conv_DG,'x-k','DisplayName','conv'); hold on;
+
+title('$\nu=0.01$');
+xlabel('Drive');
+ylabel('Transport');
+
+%% Nu = 5e-2
+KN = [ 1.5 1.6 1.7 1.8 1.9];
+Ginf_LD = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+conv_LD = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+Ginf_SG = [0.0e+0 0.0e+0 0.0e+0 2.0e-1 0.0e+0];
+conv_SG = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+Ginf_DG = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+figure
+semilogy(KN,Ginf_LD,'d--g','DisplayName','Landau'); hold on;
+semilogy(KN,Ginf_SG,'o-b','DisplayName','Sugama'); hold on;
+semilogy(KN,Ginf_DG,'^-r','DisplayName','Dougherty'); hold on;
+semilogy(KN,conv_LD,'x--k','DisplayName','changing params'); hold on;
+semilogy(KN,conv_SG,'x--k','DisplayName','changing params'); hold on;
+title('$\nu=0.05$');
+xlabel('Drive');
+ylabel('Transport');
+
+
+%% Nu = 1e-1
+KN = [ 1.5 1.6 1.7 1.8 1.9];
+Ginf_LD = [2.5e-2 2.5e-1 6.3e-1 1.3e+0 2.3e+0];
+conv_LD = [0.0e+0 2.5e-1 0.0e+0 1.2e+0 0.0e+0];
+Ginf_SG = [0.0e-9 1.3e-2 9.1e-2 3.0e-1 7.8e-1];
+Ginf_DG = [2.0e-4 6.4e-3 3.6e-2 3.1e-1 0.0e+0];
+
+figure
+semilogy(KN,Ginf_LD,'d--g','DisplayName','Landau'); hold on;
+semilogy(KN,Ginf_SG,'o-b','DisplayName','Sugama'); hold on;
+semilogy(KN,Ginf_DG,'^-r','DisplayName','Dougherty'); hold on;
+semilogy(KN,conv_LD,'x--k','DisplayName','changing params'); hold on;
+title('$\nu=0.1$');
+xlabel('Drive');
+ylabel('Transport');
+
+
+%% Nu = 5e-1
+KN = [ 1.5 1.6 1.7 1.8 1.9];
+Ginf_LD = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+conv_LD = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+Ginf_SG = [0.0e-9 3.9e-2 2.7e-1 6.6e-1 1.6e+0];
+conv_SG = [0.0e-9 0.0e+0 3.0e-1 0.0e+0 1.6e+0];
+Ginf_DG = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+
+figure
+semilogy(KN,Ginf_LD,'d--g','DisplayName','Landau'); hold on;
+semilogy(KN,Ginf_SG,'o-b','DisplayName','Sugama'); hold on;
+semilogy(KN,Ginf_DG,'^-r','DisplayName','Dougherty'); hold on;
+semilogy(KN,conv_LD,'x--k','DisplayName','changing params'); hold on;
+semilogy(KN,conv_SG,'x--k','DisplayName','changing params'); hold on;
+title('$\nu=0.5$');
+xlabel('Drive');
+ylabel('Transport');
+
+%% Primary mode stability region (confirmed with GENE)
+
+stable = [1.5 0.1; 1.5 0.05];
+unstable = [1.5 0.01];
+
+figure
+semilogy(stable(:,1),stable(:,2),'og'); hold on;
+semilogy(unstable(:,1),unstable(:,2),'or');
diff --git a/wk/Hallenbert.m b/wk/Hallenbert.m
new file mode 100644
index 0000000000000000000000000000000000000000..768341d0643953169f071d2e2e0d5ab051805ac7
--- /dev/null
+++ b/wk/Hallenbert.m
@@ -0,0 +1,79 @@
+addpath(genpath('../matlab')) % ... add
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PHYSICAL PARAMETERS
+NU = 0.1; % Collision frequency
+K_N = 1.5; % Density gradient drive
+K_T = 0.25*K_N; % Temperature '''
+K_E = 0.0; % Electrostat gradient
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+MU_X = 0.0; % X hyperdiffusivity
+MU_Y = 0.0; % Y ''
+KIN_E = 1; % Kinetic (1) or adiabatic (2) electron model
+%% GRID PARAMETERS
+NX = 150; % Spatial radial resolution ( = 2x radial modes)
+LX = 120; % Radial window size
+NY = 50; % Spatial azimuthal resolution (= azim modes)
+LY = 120; % Azimuthal window size
+NZ = 1; % number of perpendicular planes (parallel grid)
+P = 4;
+J = 2;
+%% GEOMETRY PARAMETERS
+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
+SG = 0; % Staggered z grids option
+%% TIME PARAMETERS
+TMAX = 2000; % Maximal time unit
+DT = 1e-3; % Time step
+SPS0D = 1; % Sampling per time unit for profiler
+SPS2D = 1; % Sampling per time unit for 2D arrays
+SPS3D = 1; % Sampling per time unit for 3D arrays
+SPS5D = 1/100; % Sampling per time unit for 5D arrays
+SPSCP = 0; % Sampling per time unit for checkpoints/10
+JOB2LOAD= -1;
+%% OPTIONS AND NAMING
+% Collision operator
+% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO = 'LD';
+GKCO = 1; % gyrokinetic operator
+ABCO = 1; % interspecies collisions
+NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
+SIMID = 'Hallenbert_nu_1e-01'; % Name of the simulation
+% SIMID = 'debug'; % Name of the simulation
+LINEARITY = 'nonlinear'; % (nonlinear, semilinear, linear)
+% INIT options
+INIT_PHI = 1; % Start simulation with a noisy phi (0= noisy moments 00)
+INIT_ZF = 0; ZF_AMP = 0.0;
+INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 'donothing';
+%% 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
+PMAXE = P; % Highest electron Hermite polynomial degree
+JMAXE = J; % Highest '' Laguerre ''
+PMAXI = P; % Highest ion Hermite polynomial degree
+JMAXI = J; % Highest '' Laguerre ''
+KERN = 0; % Kernel model (0 : GK)
+KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
+KPAR = 0.0; % Parellel wave vector component
+LAMBDAD = 0.0;
+NOISE0 = 1.0e-4;
+BCKGD0 = 0; % Init background
+TAU = 1.0; % e/i temperature ratio
+MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
+%% Setup and file management
+setup
+system('rm fort*.90');
+outfile = [BASIC.RESDIR,'out.txt'];
+disp(outfile);
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 3978302c2f38791881c0680287c75454758adcec..941a4ace3d2972ed825b7350413e982efe8e47c5 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -8,11 +8,85 @@ outfile ='';
outfile ='';
outfile ='';
outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
-% outfile ='nonlin_FCGK/nadiab_op_continue';
-% outfile ='Cyclone/150x150x30_5x3_Lx_200_Ly_150_q0_1.4_e_0.18_kN_2.22_kT_6.9_nu_5e-02_SGGK_adiabe';
+%% nu = 5e-1
+% Sugama
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_5e-01_SGGK';% also in 7x4
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_5e-01_SGGK';
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_5e-01_SGGK';% also in 7x4
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_5e-01_SGGK';
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_5e-01_SGGK';%also in 7x4
+
+%% nu = 1e-1
+% Landau
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/150x50_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/150x50_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-01_LDGK';
+
+% Sugama
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_SGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_SGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_SGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-01_SGGK';
+
+% Dougherty
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_DGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_DGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_DGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_DGGK';
+
+%% nu = 5e-2
+% outfile ='Hallenbert_nu_5e-02/200x32_11x6_L_120_kN_1.8_kT_0.45_nu_5e-02_SGGK';%For GENE benchmark % to analyse (added HD)
+% outfile ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGGK';
+
+% testing various NL closures
+% outfile ='Hallenbert_nu_5e-02/200x32_7x4_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGGK';
+outfile ='Hallenbert_nu_5e-02/200x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+
+%% nu = 1e-2
+% Landau
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-02_LDGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_LDGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_LDGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_LDGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-02_LDGK';
+
+% Sugama
+% outfile ='kobayashi_2015_fig1/150x150_5x3_L_100_kN_1.4_nu_5e-03_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.5_kT_0.375_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.6_kT_0.4_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/300x64_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_11x6_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_SGGK'; % To analyse (added HD)
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_11x6_L_120_kN_1.9_kT_0.475_nu_1e-02_SGGK';
+% Dougherty
+% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.5_kT_0.375_nu_1e-02_DGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_DGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_8x5_L_120_kN_1.6_kT_0.4_nu_0e+00_DGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_DGGK';
+
+%% nu = 5e-3
+% outfile ='Hallenbert_nu_5e-03/200x32_5x3_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_muxy_5e-2';
+% outfile ='Hallenbert_nu_5e-03/200x32_5x3_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_mux_5e-2_muy_6e-1';
+% outfile ='Hallenbert_nu_5e-03/200x32_11x6_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_mux_5e-2_muy_6e-1';
+% outfile ='Hallenbert_nu_5e-03/200x32_11x6_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_muxy_5e-2';
+
+%% nu = 0
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120x240_kN_1.6_kT_0.4_nu_0_mu_0.05';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_0_mu_0.05';
+% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_120_Ly_60_kN_2.5_eta_0.25_nu_0_muxy_5e-2';
+% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_120_Ly_60_kN_2.5_eta_0.25_nu_0_muxy_1e-1';
+% outfile ='Hallenbert_fig2b/200x32_11x6_Lx_120_Ly_60_kN_2.5_eta_0.25_nu_0_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nu_0_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nu_0_muxy_5e-2';
+
LOCALDIR = ['../results/',outfile,'/'];
MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
system(['mkdir -p ',MISCDIR]);
@@ -23,15 +97,15 @@ outfile ='';
outfile ='';
outfile ='';
outfile ='';
-outfile ='';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A_new/300x300_5x3_L_120_kN_1.6667_nu_1e-01_SGGK/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt';
% BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
end
%% Load the results
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 20;
+JOBNUMMIN = 00; JOBNUMMAX = 10;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
@@ -42,22 +116,33 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 11000; TAVG_1 = 11600; % Averaging times duration
-fig = plot_radial_transport_and_shear(data,TAVG_0,TAVG_1);
+TAVG_0 =100; TAVG_1 = 80000; % Averaging times duration
+% chose your field to plot in spacetime diag (uzf,szf,Gx)
+fig = plot_radial_transport_and_spacetime(data,TAVG_0,TAVG_1,'phi');
save_figure(data,fig)
end
+if 0
+%% ZF caracteristics (space time diagrams
+TAVG_0 = 500; TAVG_1 = 18000; % Averaging times duration
+% chose your field to plot in spacetime diag (uzf,szf,Gx)
+fig = ZF_spacetime(data,TAVG_0,TAVG_1);
+save_figure(data,fig)
+end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
-options.NAME = '\phi^{NZ}';
+options.NAME = '\phi';
+% options.NAME = 'v_y';
+% options.NAME = 'n_i^{NZ}';
+% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'kxky';
+options.PLAN = 'xy';
options.COMP = 1;
-options.TIME = 1500:10:5500;
+options.TIME =200:3:800;
% options.TIME = 140:0.5:160;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
@@ -68,13 +153,15 @@ if 0
% Options
options.INTERP = 0;
options.POLARPLOT = 0;
+options.AXISEQUAL = 0;
+% options.NAME = '\phi';
% options.NAME = '\Gamma_x';
-% options.NAME = 'n_i^{NZ}|';
-options.NAME = 'n_i';
-options.PLAN = 'kxz';
-options.COMP = 1;
-options.TIME = [50 70 200 500];
-data.a = data.EPS * 1000;
+options.NAME = 'n_i^{NZ}';
+% options.NAME = 'k^2n_e';
+options.PLAN = 'kxky';
+options.COMP = 'avg';
+options.TIME = 1.2e4+[600 800 1000];
+data.a = data.EPS * 2000;
fig = photomaton(data,options);
save_figure(data,fig)
end
@@ -83,10 +170,11 @@ if 0
%% 3D plot on the geometry
options.INTERP = 1;
options.NAME = 'n_i';
-options.PLANES = 16;
-options.TIME = 50;
+options.PLANES = 1;
+options.TIME = 1.2e4+[0 500 1000];
data.rho_o_R = 1e-3; % Sound larmor radius over Machine size ratio
-FIGURE = show_geometry(data,options);
+fig = show_geometry(data,options);
+save_figure(data,fig)
end
if 0
@@ -171,7 +259,7 @@ if 0
%% Zonal profiles (ky=0)
% Chose the field to plot
-FIELD = Ne00.*conj(Ne00); FNAME = 'Ne002'; FIELDLTX = '|N_e^{00}|^2'
+FIELD = data.Ne00.*conj(data.Ne00); FNAME = 'Ne002'; FIELDLTX = '|N_e^{00}|^2'
% FIELD = Ni00.*conj(Ni00); FNAME = 'Ni002'; FIELDLTX = '|N_i^{00}|^2'
% FIELD = abs(PHI); FNAME = 'absPHI'; FIELDLTX = '|\tilde\phi|';
% FIELD = PHI.*conj(PHI); FNAME = 'PHI2'; FIELDLTX = '|\tilde\phi|^2';
@@ -189,13 +277,13 @@ TNAME = [];
fig = figure; FIGNAME = ['Zonal_',FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 600,400])
plt_2 = @(x) (fftshift(x,2));
for i_ = 1:numel(tf)
- [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))];
- DATA = plt_2(squeeze(plt(FIELD,it)));
- semilogy(kx(2:end),DATA,'-','DisplayName',sprintf('$t c_s/R=%.0f$',Ts3D(it))); hold on; grid on;
+ [~,it] = min(abs(data.Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(data.Ts3D(it))];
+ d2plot = plt_2(squeeze(plt(FIELD,it)));
+ semilogy(data.kx(2:end),d2plot,'-','DisplayName',sprintf('$t c_s/R=%.0f$',data.Ts3D(it))); hold on; grid on;
xlabel(latexize(XNAME));
end
title(['$',FIELDLTX,'$ Zonal Spectrum']); legend('show');
-save_figure
+
end
if 0
diff --git a/wk/benchmark_HeLaZ_molix.m b/wk/benchmark_HeLaZ_molix.m
index bb6da7ab40c14a323db5dcb13abca47c5e88e306..053ab242ad725a04bb26c0cefac5802f19d83a88 100644
--- a/wk/benchmark_HeLaZ_molix.m
+++ b/wk/benchmark_HeLaZ_molix.m
@@ -7,7 +7,7 @@ cd ..
system('make');
cd wk
SIMDIR = '../results/benchmarks/shearless_molix_Kin_e/';
-system(['cd ',SIMDIR,';',' ./../../../bin/helaz 0; cd ../../../wk'])
+system(['cd ',SIMDIR,';',' ./../../../bin/helaz3.04 0; cd ../../../wk'])
% Run molix
% cd ../../molix
diff --git a/wk/evolve_tracers.m b/wk/evolve_tracers.m
new file mode 100644
index 0000000000000000000000000000000000000000..05e3613d6657849f92b3987d30bad59c39824a98
--- /dev/null
+++ b/wk/evolve_tracers.m
@@ -0,0 +1,260 @@
+% Options
+SHOW_FILM = 0;
+U_TIME =11000; % >0 for frozen velocity at a given time, -1 for evolving field
+Evolve_U = 1; % 0 for frozen velocity at a given time, 1 for evolving field
+Tfin = 1000;
+dt_ = 0.1;
+Nstep = ceil(Tfin/dt_);
+% Init tracers
+Np = 40; %number of tracers
+% color = tcolors;
+color = jet(Np);
+tcolors = distinguishable_colors(Np); %Their colors
+
+Na = 10000; %length of trace
+
+Traj_x = zeros(Np,Nstep);
+Traj_y = zeros(Np,Nstep);
+Disp_x = zeros(Np,Nstep);
+Disp_y = zeros(Np,Nstep);
+
+xmax = max(data.x); xmin = min(data.x);
+ymax = max(data.y); ymin = min(data.y);
+
+% Xp = 0.8*xmax*(0.5-rand(Np));
+% Yp = 0.8*ymax*(0.5-rand(Np));
+dp_ = (xmax-xmin)/(Np-1);
+Xp = linspace(xmin+dp_/2,xmax-dp_/2,Np);
+Yp = zeros(1,Np);
+Zp = zeros(Np);
+
+% position grid and velocity field
+[YY_, XX_ ,ZZ_] = meshgrid(data.y,data.x,data.z);
+[KY,KX] = meshgrid(data.ky,data.kx);
+Ux = zeros(size(XX_));
+Uy = zeros(size(XX_));
+Uz = zeros(size(XX_));
+ni = zeros(size(XX_));
+
+[~,itu_] = min(abs(U_TIME-data.Ts3D));
+% computing the frozen velocity field
+for iz = 1:data.Nz
+ Ux(:,:,iz) = real(ifft2( 1i*KY.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
+ Uy(:,:,iz) = real(ifft2(-1i*KX.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
+end
+
+
+%
+%Time loop
+t_ = 0;
+it = 1;
+itu_old = 0;
+nbytes = fprintf(2,'frame %d/%d',it,Nstep);
+if SHOW_FILM
+ fig = figure;
+end
+while(t_<Tfin)
+ if Evolve_U
+ [~,itu_] = min(abs(U_TIME+t_-data.Ts3D));
+ end
+ if Evolve_U && (itu_old ~= itu_)
+ % updating the velocity field and density field
+ for iz = 1:data.Nz
+ Ux(:,:,iz) = real(ifft2( 1i*KY.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
+ Uy(:,:,iz) = real(ifft2(-1i*KX.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
+ ni(:,:,iz) = real(ifft2(data.DENS_I(:,:,iz,itu_),data.Nx,data.Ny));
+ end
+ end
+ % evolve each tracer
+ for ip = 1:Np
+ % locate the tracer
+ % find corners of the cell
+ x_ = Xp(ip);
+ [e_x,ixC] = min(abs(x_-data.x));
+ if e_x == 0 % on the face
+ ix0 = ixC;
+ ix1 = ixC;
+ elseif x_ > data.x(ixC) % right from grid point
+ ix0 = ixC;
+ ix1 = ixC+1;
+ else % left
+ ix0 = ixC-1;
+ ix1 = ixC;
+ end
+ y_ = Yp(ip);
+ [e_y,iyC] = min(abs(y_-data.y));
+ if e_y == 0 % on the face
+ iy0 = iyC;
+ iy1 = iyC;
+ elseif y_ > data.y(iyC) % above
+ iy0 = iyC;
+ iy1 = iyC+1;
+ else % under
+ iy0 = iyC-1;
+ iy1 = iyC;
+ end
+ z_ = Zp(ip,1);
+ [e_z,izC] = min(abs(z_-data.z));
+ if e_z == 0 % on the face
+ iz0 = izC;
+ iz1 = izC;
+ elseif z_ > data.z(izC) % before
+ iz0 = izC;
+ iz1 = izC+1;
+ else % behind
+ iz0 = izC-1;
+ iz1 = izC;
+ end
+ x0 = data.x(ix0); x1 = data.x(ix1); %left right
+ y0 = data.x(iy0); y1 = data.y(iy1); %down top
+ z0 = data.z(iz0); z1 = data.z(iz1); %back front
+ if(e_x > 0)
+ ai__ = (x_ - x0)/(x1-x0); % interp coeff x
+ else
+ ai__ = 0;
+ end
+ if(e_y > 0)
+ a_i_ = (y_ - y0)/(y0-y1); % interp coeff y
+ else
+ a_i_ = 0;
+ end
+ if(e_z > 0)
+ a__i = (z_ - z0)/(z1-z0); % interp coeff z
+ else
+ a__i = 0;
+ end
+ % interp velocity and density
+ %velocity values
+ u000 = [Ux(ix0,iy0,iz0) Uy(ix0,iy0,iz0) ni(ix0,iy0,iz0)];
+ u001 = [Ux(ix0,iy0,iz1) Uy(ix0,iy0,iz1) ni(ix0,iy0,iz1)];
+ u010 = [Ux(ix0,iy1,iz0) Uy(ix0,iy1,iz0) ni(ix0,iy1,iz0)];
+ u011 = [Ux(ix0,iy1,iz1) Uy(ix0,iy1,iz1) ni(ix0,iy1,iz1)];
+ u100 = [Ux(ix1,iy0,iz0) Uy(ix1,iy0,iz0) ni(ix1,iy0,iz0)];
+ u101 = [Ux(ix1,iy0,iz1) Uy(ix1,iy0,iz1) ni(ix1,iy0,iz1)];
+ u110 = [Ux(ix1,iy1,iz0) Uy(ix1,iy1,iz0) ni(ix1,iy1,iz0)];
+ u111 = [Ux(ix1,iy1,iz1) Uy(ix1,iy1,iz1) ni(ix1,iy1,iz1)];
+ %linear interpolation
+ linterp = @(x1,x2,a) (1-a)*x1 + a*x2;
+
+ u_00 = linterp(u000,u100,ai__);
+ u_10 = linterp(u010,u110,ai__);
+ u__0 = linterp(u_00,u_10,a_i_);
+
+ u_01 = linterp(u001,u101,ai__);
+ u_11 = linterp(u011,u111,ai__);
+ u__1 = linterp(u_01,u_11,a_i_);
+
+ u___ = linterp(u__0,u__1,a__i);
+
+ % push the particle
+ q = sign(-u___(3));
+% q =1;
+ x_ = x_ + dt_*u___(1)*q;
+ y_ = y_ + dt_*u___(2)*q;
+
+ % apply periodic boundary conditions
+ if(x_ > xmax)
+ x_ = xmin + (x_ - xmax);
+ elseif(x_ < xmin)
+ x_ = xmax + (x_ - xmin);
+ end
+ if(y_ > ymax)
+ y_ = ymin + (y_ - ymax);
+ elseif(y_ < ymin)
+ y_ = ymax + (y_ - ymin);
+ end
+ % store data
+ % Trajectory
+ Traj_x(ip,it) = x_;
+ Traj_y(ip,it) = y_;
+ % Displacement
+ Disp_x(ip,it) = Disp_x(ip,it) + dt_*u___(1)*q;
+ Disp_y(ip,it) = Disp_y(ip,it) + dt_*u___(2)*q;
+ % update position
+ Xp(ip) = x_; Yp(ip) = y_;
+ end
+ %% Movie
+ if Evolve_U && (itu_old ~= itu_) && SHOW_FILM
+ % updating the velocity field
+ clf(fig);
+ F2P = real(ifft2(data.PHI(:,:,iz,itu_),data.Nx,data.Ny));
+ scale = max(max(abs(F2P))); % Scaling to normalize
+ pclr = pcolor(XX_,YY_,F2P/scale);
+ colormap(bluewhitered);
+ set(pclr, 'edgecolor','none'); hold on; caxis([-5,5]);
+ for ip = 1:Np
+ ia0 = max(1,it-Na);
+ plot(Traj_x(ip,ia0:it),Traj_y(ip,ia0:it),'.','Color',color(ip,:)); hold on
+ end
+ for ip = 1:Np
+ plot(Traj_x(ip, 1),Traj_y(ip, 1),'xk'); hold on
+ plot(Traj_x(ip,it),Traj_y(ip,it),'ok','MarkerFaceColor',color(ip,:)); hold on
+ end
+ title(['$t \approx$', sprintf('%.3d',ceil(data.Ts3D(itu_)))]);
+ axis equal
+ xlim([xmin xmax]); ylim([ymin ymax]);
+ drawnow
+ end
+ t_ = t_ + dt_; it = it + 1; itu_old = itu_;
+ % terminal info
+ while nbytes > 0
+ fprintf('\b')
+ nbytes = nbytes - 1;
+ end
+ nbytes = fprintf(2,'frame %d/%d',it,Nstep);
+end
+Nt = it;
+%% Plot trajectories and statistics
+xtot = Disp_x;
+ytot = Disp_y;
+% aggregate the displacements
+for iq = 1:Np
+ x_ = 0; y_ = 0;
+ for iu = 1:Nt-1
+ x_ = x_ + Disp_x(iq,iu);
+ y_ = y_ + Disp_y(iq,iu);
+ xtot(iq,iu) = x_;
+ ytot(iq,iu) = y_;
+ end
+end
+ma = @(x) x;%movmean(x,100);
+time_ = linspace(U_TIME,U_TIME+Tfin,it-1);
+figure;
+subplot(221);
+for ip = 1:Np
+% plot(ma(time_),ma(Traj_x(ip,:)),'-','Color',tcolors(ip,:)); hold on
+ plot(ma(time_),xtot(ip,:),'-','Color',tcolors(ip,:)); hold on
+% plot(ma(time_),ma(Disp_x(ip,:)),'-','Color',tcolors(ip,:)); hold on
+end
+ylabel('$x_p$');
+xlim(U_TIME + [0 Tfin]);
+
+subplot(222);
+ plot(time_,mean(xtot,1)); hold on
+ fit = polyfit(time_,mean(xtot,1),1);
+ plot(time_,fit(1)*time_+fit(2),'--k','DisplayName',['$\alpha=',num2str(fit(1)),'$']); hold on
+ylabel('$\langle x \rangle_p$');
+xlim(U_TIME + [0 Tfin]);
+
+subplot(223);
+for ip = 1:Np
+% plot(ma(time_),ma(Traj_y(ip,:)),'-','Color',tcolors(ip,:)); hold on
+ plot(ma(time_),ytot(ip,:),'-','Color',tcolors(ip,:)); hold on
+% plot(ma(time_),ma(Disp_y(ip,:)),'-','Color',tcolors(ip,:)); hold on
+end
+xlabel('time');
+ylabel('$y_p$');
+xlim(U_TIME + [0 Tfin]);
+
+subplot(224);
+ plot(time_,mean(ytot,1)); hold on
+xlabel('time');
+ylabel('$\langle y \rangle_p$');
+xlim(U_TIME + [0 Tfin]);
+
+if 0
+ %%
+ figure
+% hist(reshape(xtot,[Np*(Nt-1) 1],Np))
+ hist(xtot(:,2000),Np/20)
+end
diff --git a/wk/linear_1D_entropy_mode.m b/wk/linear_1D_entropy_mode.m
index 5c23f2f39bcc3b54a50ba30b8d70dc96be6876df..a6dd80cb252bfe9e60212251b99ee403c3bae014 100644
--- a/wk/linear_1D_entropy_mode.m
+++ b/wk/linear_1D_entropy_mode.m
@@ -8,14 +8,14 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_N = 1/0.6; % Density gradient drive
-K_T = 0.0; % Temperature '''
+K_N = 1.8; % Density gradient drive
+K_T = 0.25*K_N; % Temperature '''
K_E = 0.0; % Electrostat '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
%% GRID PARAMETERS
-NX = 150; % real space x-gridpoints
+NX = 64; % real space x-gridpoints
NY = 1; % '' y-gridpoints
-LX = 200; % Size of the squared frequency domain
+LX = 100; % Size of the squared frequency domain
LY = 1; % Size of the squared frequency domain
NZ = 1; % number of perpendicular planes (parallel grid)
Q0 = 1.0; % safety factor
@@ -23,8 +23,8 @@ SHEAR = 0.0; % magnetic shear
EPS = 0.0; % inverse aspect ratio
SG = 1; % Staggered z grids option
%% TIME PARMETERS
-TMAX = 200; % Maximal time unit
-DT = 1e-2; % Time step
+TMAX = 100; % Maximal time unit
+DT = 1e-2c; % 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
@@ -32,19 +32,19 @@ SPS5D = 1; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
-SIMID = 'Linear_entropy_mode'; % Name of the simulation
+SIMID = 'Linear_Hallenbert'; % Name of the simulation
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
KIN_E = 1;
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO = 'LD';
-GKCO = 1; % gyrokinetic operator
+CO = 'SG';
+GKCO = 0; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
-INIT_PHI= 0; % Start simulation with a noisy phi
+INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
@@ -56,11 +56,13 @@ W_NAPJ = 1; W_SAPJ = 0;
% unused
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
kmax = NX*pi/LX;% Highest fourier mode
-NU_HYP = 0.0; % Hyperdiffusivity coefficient
-MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
+MU = 0.05; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
-MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
+MU_X = MU; %
+MU_Y = MU; %
+MU_Z = 0.0; %
+MU_P = 0.0; %
+MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
@@ -70,10 +72,10 @@ CURVB = 1.0;
if 1
%% Parameter scan over PJ
-% PA = [2 4];
-% JA = [1 2];
-PA = [4];
-JA = [2];
+% PA = [20];
+% JA = [10];
+PA = [4 6 8];
+JA = [2 3 4];
DTA= DT*ones(size(JA));%./sqrt(JA);
% DTA= DT;
mup_ = MU_P;
@@ -118,7 +120,7 @@ for i = 1:Nparam
if 0
%% Fit verification
figure;
- for i = 1:1:Nx/2+1
+ for i = 1:1:NX/2+1
X_ = Ts3D(:); Y_ = squeeze(abs(Ni00(i,1,1,:)));
semilogy(X_,Y_,'DisplayName',['k=',num2str(kx(i))]); hold on;
end
@@ -131,7 +133,9 @@ fig = figure; FIGNAME = 'linear_study';
plt = @(x) x;
% subplot(211)
for i = 1:Nparam
- clr = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:);
+% colors = line_colors;
+ colors = jet(Nparam);
+ clr = colors(mod(i-1,numel(line_colors(:,1)))+1,:);
linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1);
plot(plt(SCALE*kx),plt(gamma_phi(i,1:end)),...
'Color',clr,...
diff --git a/wk/linear_1D_damping.m b/wk/linear_damping.m
similarity index 77%
rename from wk/linear_1D_damping.m
rename to wk/linear_damping.m
index 35e2a5313cf702586c0cf1280e67fd949c11492b..394e93b56e88d0c57025e63a28d3d263b247502c 100644
--- a/wk/linear_1D_damping.m
+++ b/wk/linear_damping.m
@@ -8,40 +8,43 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
TAU = 1.0; % e/i temperature ratio
+K_N = 0.0; % Density gradient drive
+K_T = 0.0; % Temperature '''
+K_E = 0.0; % Electrostat '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
%% GRID PARAMETERS
-NX = 150; % real space x-gridpoints
+NX = 64; % real space x-gridpoints
NY = 1; % '' y-gridpoints
-LX = 200; % Size of the squared frequency domain
+LX = 100; % Size of the squared frequency domain
LY = 1; % Size of the squared frequency domain
NZ = 1; % number of perpendicular planes (parallel grid)
Q0 = 1.0; % safety factor
SHEAR = 0.0; % magnetic shear
EPS = 0.0; % inverse aspect ratio
-SG = 0; % Staggered z grids option
+SG = 1; % Staggered z grids option
%% TIME PARMETERS
-TMAX = 20; % Maximal time unit
-DT = 1e-2; % Time step
+TMAX = 500; % Maximal time unit
+DT = 2e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
-SPS3D = 5; % Sampling per time unit for 2D arrays
-SPS5D = 1/5; % Sampling per time unit for 5D arrays
+SPS3D = 1; % Sampling per time unit for 2D arrays
+SPS5D = 1; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
-NOISE0 = 0.0; % Init noise amplitude
-BCKGD0 = 1.0; % Init background
-SIMID = 'Linear_damping'; % Name of the simulation
-LINEARITY = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
+SIMID = 'damping_analysis'; % Name of the simulation
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
KIN_E = 1;
% Collision operator
-% (0:L.Bernstein, 1:Dougherty, 2:Sugama, 3:Pitch angle, 4:Full Couloumb ; +/- for GK/DK)
-CO = 4;
+% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO = 'DG';
+GKCO = 1; % gyrokinetic operator
+ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
-INIT_PHI= 0; % Start simulation with a noisy phi
+INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
@@ -56,22 +59,24 @@ kmax = NX*pi/LX;% Highest fourier mode
NU_HYP = 0.0; % Hyperdiffusivity coefficient
MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
-MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
+MU_X = 0.0; %
+MU_Y = 0.0; %
+MU_Z = 0.0; %
+MU_P = 0.0; %
+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_N = 0.0; % Density gradient drive
-K_T = 0.0; % Temperature '''
-K_E = 0.0; % Electrostat '''
%% PARAMETER SCANS
if 1
%% Parameter scan over PJ
-% PA = [2 4];
-% JA = [1 2];
-PA = [2];
-JA = [1];
+PA = [4];
+JA = [2];
+% PA = [2 4 8];
+% JA = [1 2 4];
DTA= DT*ones(size(JA));%./sqrt(JA);
% DTA= DT;
mup_ = MU_P;
@@ -91,34 +96,36 @@ for i = 1:Nparam
% Run linear simulation
if RUN
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz3 1 6 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ./../../../bin/helaz3 1 2 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz3 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ./../../../bin/helaz 1 2 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
end
% Load and process results
%%
filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5'];
load_results
+ tfit = []; gfit = [];
for ikx = 1:NX/2+1
- tend = 2;%max(Ts3D(abs(Ni00(ikx,1,1,:))~=0));
- tstart = 0;
+ tend = max(Ts3D);
+ tstart = 0.5*tend;
[~,itstart] = min(abs(Ts3D-tstart));
[~,itend] = min(abs(Ts3D-tend));
trange = itstart:itend;
% exp fit on moment 00
X_ = Ts3D(trange); Y_ = squeeze(abs(Ni00(ikx,1,1,trange)));
- gamma_Ni00(i,ikx) = LinearFit_s(X_,Y_);
+ [gamma_Ni00(i,ikx),tfit(ikx,:),gfit(ikx,:)] = LinearFit_s(X_,Y_);
% exp fit on phi
X_ = Ts3D(trange); Y_ = squeeze(abs(PHI(ikx,1,1,trange)));
gamma_phi (i,ikx) = LinearFit_s(X_,Y_);
end
gamma_Ni00(i,:) = real(gamma_Ni00(i,:));% .* (gamma_Ni00(i,:)>=0.0));
gamma_Nipj(i,:) = real(gamma_Nipj(i,:));% .* (gamma_Nipj(i,:)>=0.0));
- if 0
+ if 1
%% Fit verification
figure;
- for i = 1:1:NX/2+1
- X_ = Ts3D(:); Y_ = squeeze(abs(Ni00(i,1,1,:)));
- semilogy(X_,Y_,'DisplayName',['k=',num2str(kx(i))]); hold on;
+ for i = 1:5:NX/2+1
+% X_ = Ts3D(:); Y_ = squeeze(abs(Ni00(i,1,1,:)));
+% semilogy(X_,Y_,'DisplayName',['k=',num2str(kx(i))]); hold on;
+ plot(tfit(ikx,:),gfit(i,:),'DisplayName',['k=',num2str(kx(i))]); hold on;
end
end
diff --git a/wk/local_run.m b/wk/local_run.m
index c5d950fe21d2ed88de46e7f4262a6d8b123874ea..f2581ae23c035f6a7522a169cb2bc6e7dc187f2e 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -4,18 +4,18 @@ addpath(genpath('../matlab')) % ... add
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
-K_N = 1/0.6; % Density gradient drive
+NU = 0.005; % Collision frequency
+K_N = 1.4; % Density gradient drive
K_T = 0.0; % Temperature '''
K_E = 0.0; % Electrostat gradient
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
NU_HYP = 0.0;
KIN_E = 1; % Kinetic (1) or adiabatic (2) electron model
%% GRID PARAMETERS
-NX = 200; % Spatial radial resolution ( = 2x radial modes)
-LX = 120; % Radial window size
-NY = 100; % Spatial azimuthal resolution (= azim modes)
-LY = 120; % Azimuthal window size
+NX = 150; % Spatial radial resolution ( = 2x radial modes)
+LX = 100; % Radial window size
+NY = 150; % Spatial azimuthal resolution (= azim modes)
+LY = 100; % Azimuthal window size
NZ = 1; % number of perpendicular planes (parallel grid)
P = 4;
J = 2;
@@ -27,24 +27,26 @@ GRADB = 1.0; % Magnetic gradient
CURVB = 1.0; % Magnetic curvature
SG = 0; % Staggered z grids option
%% TIME PARAMETERS
-TMAX = 1000; % Maximal time unit
-DT = 1e-2; % Time step
+TMAX = 120; % Maximal time unit
+DT = 2e-2; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS3D = 1/2; % Sampling per time unit for 3D arrays
+SPS3D = 1; % Sampling per time unit for 3D arrays
SPS5D = 1/200; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints/10
-JOB2LOAD= 1;
+JOB2LOAD= -1;
%% OPTIONS AND NAMING
% Collision operator
-% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; 4 : Coulomb; +/- for GK/DK)
-CO = 2;
-CLOS = 0; % Closure model (0: =0 truncation)
+% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO = 'SG';
+GKCO = 1; % gyrokinetic operator
+ABCO = 1; % interspecies collisions
NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-SIMID = 'semi_linear'; % Name of the simulation
-LINEARITY = 'nonlinear'; % (nonlinear, semilinear, linear)
+SIMID = 'kobayashi_2015_fig1'; % Name of the simulation
+% SIMID = 'debug'; % Name of the simulation
+LINEARITY = 'linear'; % (nonlinear, semilinear, linear)
% INIT options
-INIT_PHI = 0; % Start simulation with a noisy phi (0= noisy moments 00)
+INIT_PHI = 1; % Start simulation with a noisy phi (0= noisy moments 00)
INIT_ZF = 0; ZF_AMP = 0.0;
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 'donothing';
%% OUTPUTS
@@ -67,8 +69,8 @@ LAMBDAD = 0.0;
kmax = NX*pi/LX;% Highest fourier mode
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
-NOISE0 = 1.0e-5;
-BCKGD0 = 0.0; % Init background
+NOISE0 = 0;%1.0e-4;
+BCKGD0 = 1e-4; % Init background
TAU = 1.0; % e/i temperature ratio
MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index f9836ed232d0504c7f18113acc5009525ff34986..841237cb613151baef69e6fb9b1e29cc66dd7f32 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -3,7 +3,8 @@ addpath(genpath('../matlab')) % ... add
SUBMIT = 1; % To submit the job automatically
CHAIN = 2; % To chain jobs (CHAIN = n will launch n jobs in chain)
% EXECNAME = 'helaz3_dbg';
- EXECNAME = 'helaz3';
+ EXECNAME = 'helaz3.03';
+ SIMID = 'simulation_A_new';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -35,25 +36,24 @@ P = 4;
J = 2;
%% TIME PARAMETERS
TMAX = 10000; % Maximal time unit
-DT = 2e-3; % Time step
+DT = 7e-3; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
SPS3D = 1/2; % Sampling per time unit for 3D arrays
SPS5D = 1/50; % Sampling per time unit for 5D arrays
-JOB2LOAD= -1; % start from t=0 if <0, else restart from outputs_$job2load
+JOB2LOAD= 0; % start from t=0 if <0, else restart from outputs_$job2load
%% OPTIONS AND NAMING
% Collision operator
-% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK)
-CO = 2;
+% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO = 'LR';
+GKCO = 1; % gyrokinetic operator
+ABCO = 1; % interspecies collisions
CLOS = 0; % Closure model (0: =0 truncation)
-NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-% SIMID = 'test_chained_job'; % Name of the simulation
-SIMID = 'simulation_A'; % Name of the simulation
-% SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
-LINEARITY = 1; % activate non-linearity (is cancelled if KXEQ0 = 1)
+NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
+LINEARITY = 'nonlinear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% INIT options
INIT_ZF = 0; ZF_AMP = 0.0;
-INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
+INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 'donothing';
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
@@ -111,6 +111,8 @@ if(SUBMIT)
SBATCH_CMD = ['sbatch --dependency=afterok:',jobid_,' batch_script.sh\n'];
disp(SBATCH_CMD);
JOB2LOAD= JOB2LOAD+1;
+% DT = 1e-1;
+% NL_CLOS = -1;
setup
write_sbash_marconi
[~,job_info_] = system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh');
@@ -118,5 +120,5 @@ if(SUBMIT)
end
end
end
-system('rm fort*.90');
+% system('rm fort*.90');
disp('done');
diff --git a/wk/mode_growth_meter.m b/wk/mode_growth_meter.m
index e00a08aab3eb88eda7a56665088929357f3993bb..3b90ca720d00b0c61e4217d0268a39899b29f78e 100644
--- a/wk/mode_growth_meter.m
+++ b/wk/mode_growth_meter.m
@@ -1,12 +1,18 @@
-MODES_SELECTOR = 1; %(1:Zonal, 2: NZonal, 3: ky=kx)
-NORMALIZED = 1;
-options.K2PLOT = 0.01:0.05:1.0;
-options.TIME =500:0.5:1000;
+MODES_SELECTOR = 2; %(1:Zonal, 2: NZonal, 3: ky=kx)
+NORMALIZED = 0;
+options.K2PLOT = 0.01:0.01:1.0;
+options.TIME =4000:1:4200;
DATA = data;
OPTIONS = options;
t = OPTIONS.TIME;
iz = 1;
+[~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(:,1,1,:))),2)));
+
+FRAMES = zeros(size(OPTIONS.TIME));
+for i = 1:numel(OPTIONS.TIME)
+ [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D));
+end
if MODES_SELECTOR == 1
if NORMALIZED
@@ -37,16 +43,11 @@ elseif MODES_SELECTOR == 3
MODESTR = 'Diag modes';
end
-
-FRAMES = zeros(size(OPTIONS.TIME));
-for i = 1:numel(OPTIONS.TIME)
- [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D));
-end
-
-MODES = zeros(size(OPTIONS.K2PLOT));
-for i = 1:numel(OPTIONS.K2PLOT)
- [~,MODES(i)] =min(abs(OPTIONS.K2PLOT(i)-k));
-end
+MODES = 1:numel(k);
+% MODES = zeros(size(OPTIONS.K2PLOT));
+% for i = 1:numel(OPTIONS.K2PLOT)
+% [~,MODES(i)] =min(abs(OPTIONS.K2PLOT(i)-k));
+% end
% plt = @(x,ik) abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES))));
@@ -74,7 +75,10 @@ subplot(1,3,2)
mod2plot = round(linspace(2,numel(MODES),15));
for i_ = mod2plot
semilogy(t,plt(DATA.PHI,MODES(i_))); hold on;
- semilogy(t,exp(gamma(i_).*t+amp(i_)),'--k')
+% semilogy(t,exp(gamma(i_).*t+amp(i_)),'--k')
+ end
+ if MODES_SELECTOR == 1
+ semilogy(t,plt(DATA.PHI,ikzf),'--k');
end
xlim([t(1) t(end)]); %ylim([1e-5 1])
xlabel('$t c_s /\rho_s$'); ylabel(['$|\phi_{',kstr,'}|$']);
@@ -83,5 +87,8 @@ subplot(1,3,2)
subplot(1,3,3)
plot(k(MODES),gamma); hold on;
plot(k(MODES(mod2plot)),gamma(mod2plot),'x')
+ if MODES_SELECTOR == 1
+ plot(k(ikzf),gamma(ikzf),'ok');
+ end
xlabel(['$',kstr,'\rho_s$']); ylabel('$\gamma$');
title('Growth rates')
\ No newline at end of file
diff --git a/wk/quick_gene_load.m b/wk/quick_gene_load.m
new file mode 100644
index 0000000000000000000000000000000000000000..40cb7617192efc0128a863a2291c34ee8c9916cf
--- /dev/null
+++ b/wk/quick_gene_load.m
@@ -0,0 +1,36 @@
+%% gyroLES plot
+% genepath = '/misc/gene_results';
+% % simpath = '/NL_Zpinch_Kn_1.7_eta_0.25_nuSG_1e-1_gyroLES/';
+% simpath = '/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_gyroLES/';
+% filename = 'GyroLES_ions.dat';
+% toload = [genepath simpath filename];
+% data = readtable(toload);
+% t = data.Var1;
+% mu_x = data.Var2;
+% mu_y = data.Var3;
+% figure;
+% plot(t,mu_x); hold on
+% plot(t,mu_y);
+% legend('mu_x','mu_y'); xlabel('time'); ylabel('Hyper-diff');
+% title('gyroLES results for H&P fig. 2c')
+
+%% Plot linear results
+path = '/home/ahoffman/gene/linear_zpinch_results/';
+% fname ='tmp.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSG_0.047_32x16.txt';
+% fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.047_32x16.txt';
+% fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.0235_32x16.txt';
+% fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.0047_32x16.txt';
+% fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.0047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.4_KT_0.35_nuSG_0.0047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.6_KT_0.4_nuSG_0.0047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.45_nu_0_32x16.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSG_0.0235_64x32.txt';
+% fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.235_64x32.txt';
+% fname ='GENE_LIN_Kn_1.7_KT_0.425_nuSG_0.235_64x32.txt';
+fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSGDK_0.047_32x16.txt';
+data_ = load([path,fname]);
+
+figure
+plot(data_(:,2),data_(:,3),'-dk','DisplayName','GENE');
\ No newline at end of file
diff --git a/wk/tmp.txt b/wk/tmp.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391