Skip to content
Snippets Groups Projects
Commit 44d4d8b2 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

script save

parent 2bf9db7d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 423 additions and 102 deletions
matlab/load/quick_gene_load.m
+ 2
− 2
View file @ 44d4d8b2
......@@ -73,8 +73,8 @@ path = '/home/ahoffman/gene/linear_CBC_results/';
%----------Convergence nvpar shearless CBC
% fname = 'CBC_salpha_nz_24_nv_scan_nw_16_adiabe.txt';
% fname = 'CBC_miller_nz_24_nv_scan_nw_16_adiabe.txt';
% fname = 'CBC_salpha_nz_24_nv_scan_nw_16_kine.txt';
fname = 'CBC_miller_nz_24_nv_scan_nw_16_kine.txt';
fname = 'CBC_salpha_nz_24_nv_scan_nw_16_kine.txt';
% fname = 'CBC_miller_nz_24_nv_scan_nw_16_kine.txt';
%---------- CBC
% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_kine.txt';
......
matlab/plot/create_film.m
+ 1
− 1
View file @ 44d4d8b2
......@@ -12,7 +12,7 @@ switch OPTIONS.PLAN
case 'RZ'
toplot = poloidal_plot(DATA,OPTIONS);
otherwise
toplot = process_field(DATA,OPTIONS);
toplot = process_field_v2(DATA,OPTIONS);
end
%%
FILENAME = [DATA.localdir,toplot.FILENAME,format];
......
matlab/plot/photomaton.m
+ 8
− 3
View file @ 44d4d8b2
......@@ -5,7 +5,7 @@ switch OPTIONS.PLAN
case 'RZ'
toplot = poloidal_plot(DATA,OPTIONS);
otherwise
toplot = process_field(DATA,OPTIONS);
toplot = process_field_v2(DATA,OPTIONS);
end
FNAME = toplot.FILENAME;
FRAMES = toplot.FRAMES;
......@@ -15,9 +15,14 @@ Ncols = ceil(Nframes/Nrows);
%
TNAME = [];
FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows])
frame_max = max(max(max(abs(toplot.FIELD(:,:,:)))));
for i_ = 1:numel(FRAMES)
subplot(Nrows,Ncols,i_); TNAME = [TNAME,'_',sprintf('%.0f',DATA.Ts3D(FRAMES(i_)))];
scale = max(max(abs(toplot.FIELD(:,:,i_)))); % Scaling to normalize
if OPTIONS.NORMALIZE
scale = frame_max; % Scaling to normalize
else
scale = 1;
end
if ~strcmp(OPTIONS.PLAN,'sx')
tshot = DATA.Ts3D(FRAMES(i_));
pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale); set(pclr, 'edgecolor','none');
......@@ -25,7 +30,7 @@ FIGURE.fig = figure; set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows])
pbaspect(toplot.ASPECT)
end
if ~strcmp(OPTIONS.PLAN,'kxky')
caxis([-1,1]);
caxis([-1,1]*frame_max/scale);
colormap(bluewhitered);
if OPTIONS.INTERP
shading interp;
......
matlab/process_field_v2.m 0 → 100644
+ 258
− 0
View file @ 44d4d8b2
function [ TOPLOT ] = process_field( DATA, OPTIONS )
%UNTITLED4 Summary of this function goes here
% Detailed explanation goes here
%% Setup time
%%
FRAMES = zeros(size(OPTIONS.TIME));
for i = 1:numel(OPTIONS.TIME)
[~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D));
end
FRAMES = unique(FRAMES);
%% Setup the plot geometry
[KX, KY] = meshgrid(DATA.kx, DATA.ky);
directions = {'x','y','z'};
Nx = DATA.Nx; Ny = DATA.Ny; Nz = DATA.Nz; Nt = numel(FRAMES);
POLARPLOT = OPTIONS.POLARPLOT;
LTXNAME = OPTIONS.NAME;
switch OPTIONS.PLAN
case 'xy'
XNAME = '$x$'; YNAME = '$y$';
[X,Y] = meshgrid(DATA.x,DATA.y);
REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'xz'
XNAME = '$x$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.x);
REALP = 1; COMPDIM = 1; SCALE = 0;
case 'yz'
XNAME = '$y$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.y);
REALP = 1; COMPDIM = 2; SCALE = 0;
case 'kxky'
XNAME = '$k_x$'; YNAME = '$k_y$';
[X,Y] = meshgrid(DATA.kx,DATA.ky);
REALP = 0; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'kxz'
XNAME = '$k_x$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.kx);
REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0;
case 'kyz'
XNAME = '$k_y$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.ky);
REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
end
Xmax_ = max(max(abs(X))); Ymax_ = max(max(abs(Y)));
Lmin_ = min([Xmax_,Ymax_]);
Rx = Xmax_/Lmin_ * SCALE + (1-SCALE)*1.2;
Ry = Ymax_/Lmin_ * SCALE + (1-SCALE)*1.2;
ASPECT = [Rx, Ry, 1];
DIMENSIONS = [500, 1000, OPTIONS.RESOLUTION*Rx, OPTIONS.RESOLUTION*Ry];
% Polar grid
POLARNAME = [];
if POLARPLOT
POLARNAME = 'polar';
X__ = (X+DATA.a).*cos(Y);
Y__ = (X+DATA.a).*sin(Y);
X = X__;
Y = Y__;
XNAME='X';
YNAME='Z';
DIMENSIONS = [100, 100, OPTIONS.RESOLUTION, OPTIONS.RESOLUTION];
ASPECT = [1,1,1];
sz = size(X);
FIELD = zeros(sz(1),sz(2),Nt);
else
sz = size(X);
FIELD = zeros(sz(1),sz(2),Nt);
end
%% Process the field to plot
switch OPTIONS.COMP
case 'avg'
compr = @(x) mean(x,COMPDIM);
COMPNAME = ['avg',directions{COMPDIM}];
FIELDNAME = ['\langle ',LTXNAME,'\rangle_',directions{COMPDIM}];
case 'max'
compr = @(x) max(x,[],COMPDIM);
COMPNAME = ['max',directions{COMPDIM}];
FIELDNAME = ['\max_',directions{COMPDIM},' ',LTXNAME];
otherwise
switch COMPDIM
case 1
i = OPTIONS.COMP;
compr = @(x) x(i,:,:);
if REALP
COMPNAME = sprintf(['y=','%2.1f'],DATA.x(i));
else
COMPNAME = sprintf(['k_y=','%2.1f'],DATA.kx(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 2
i = OPTIONS.COMP;
compr = @(x) x(:,i,:);
if REALP
COMPNAME = sprintf(['x=','%2.1f'],DATA.y(i));
else
COMPNAME = sprintf(['k_x=','%2.1f'],DATA.ky(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 3
i = OPTIONS.COMP;
compr = @(x) x(:,:,i);
COMPNAME = sprintf(['z=','%2.1f'],DATA.z(i));
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
end
end
switch REALP
case 1 % Real space plot
INTERP = OPTIONS.INTERP;
process = @(x) real(fftshift(ifourier_GENE(x)));
shift_x = @(x) x;
shift_y = @(x) x;
case 0 % Frequencies plot
INTERP = 0;
switch COMPDIM
case 1
process = @(x) abs(fftshift(x,2));
shift_x = @(x) fftshift(x,1);
shift_y = @(x) fftshift(x,1);
case 2
process = @(x) abs((x));
shift_x = @(x) (x);
shift_y = @(x) (x);
case 3
process = @(x) abs(fftshift(x,2));
shift_x = @(x) fftshift(x,2);
shift_y = @(x) fftshift(x,2);
end
end
switch OPTIONS.NAME
case '\phi' %ES pot
NAME = 'phi';
FLD_ = DATA.PHI(:,:,:,FRAMES); % data to plot
OPE_ = 1; % Operation on data
case '\psi' %EM pot
NAME = 'psi';
FLD_ = DATA.PSI(:,:,:,FRAMES);
OPE_ = 1;
case '\phi^{NZ}' % non-zonal ES pot
NAME = 'phiNZ';
FLD_ = DATA.PHI(:,:,:,FRAMES);
OPE_ = (KY~=0);
case 'v_{Ey}' % y-comp of ExB velocity
NAME = 'vy';
FLD_ = DATA.PHI(:,:,:,FRAMES);
OPE_ = -1i*KX;
case 'v_{Ex}' % x-comp of ExB velocity
NAME = 'vx';
FLD_ = DATA.PHI(:,:,:,FRAMES);
OPE_ = -1i*KY;
case 's_{Ey}' % y-comp of ExB shear
NAME = 'sy';
FLD_ = DATA.PHI(:,:,:,FRAMES);
OPE_ = KX.^2;
case 's_{Ex}' % x-comp of ExB shear
NAME = 'sx';
FLD_ = DATA.PHI(:,:,:,FRAMES);
OPE_ = KY.^2;
case '\omega_z' % ES pot vorticity
NAME = 'vorticity';
FLD_ = DATA.PHI(:,:,:,FRAMES);
OPE_ = -(KX.^2+KY.^2);
case 'N_i^{00}' % first ion gyromoment
NAME = 'Ni00';
FLD_ = DATA.Ni00(:,:,:,FRAMES);
OPE_ = 1;
case 'N_e^{00}' % first electron gyromoment
NAME = 'Ne00';
FLD_ = DATA.Ne00(:,:,:,FRAMES);
OPE_ = 1;
case 'n_i' % ion density
NAME = 'ni';
FLD_ = DATA.DENS_I(:,:,:,FRAMES);
OPE_ = 1;
case 'n_e' % electron density
NAME = 'ne';
FLD_ = DATA.DENS_E(:,:,:,FRAMES);
OPE_ = 1;
case 'k^2n_e' % electron vorticity
NAME = 'k2ne';
FLD_ = DATA.DENS_E(:,:,:,FRAMES);
OPE_ = -(KX.^2+KY.^2);
case 'n_i-n_e' % polarisation
NAME = 'pol';
OPE_ = 1;
FLD_ = DATA.DENS_I(:,:,:,FRAMES)+DATA.DENS_E(:,:,:,FRAMES);
case 'T_i' % ion temperature
NAME = 'Ti';
FLD_ = DATA.TEMP_I(:,:,:,FRAMES);
OPE_ = 1;
case 'G_x' % ion particle flux
NAME = 'Gx';
FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
OPE_ = 1;
for it = 1:numel(FRAMES)
tmp = zeros(DATA.Ny,DATA.Nx,Nz);
for iz = 1:DATA.Nz
vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
gx_ = vx_.*ni_;
% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))),2));
tmp(:,:,iz) = abs((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))));
end
FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
end
case 'Q_x' % ion heat flux
NAME = 'Qx';
FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
OPE_ = 1;
for it = 1:numel(FRAMES)
tmp = zeros(DATA.Ny,DATA.Nx,Nz);
for iz = 1:DATA.Nz
vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
Ti_ = real((ifourier_GENE( DATA.TEMP_I(:,:,iz,FRAMES(it)))));
qx_ = vx_.*ni_.*Ti_;
% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))),2));
tmp(:,:,iz) = abs((squeeze(fft2(qx_,DATA.Ny,DATA.Nx))));
end
FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
end
otherwise
disp('Fieldname not recognized');
return
end
% Process the field according to the 2D plane and the space (real/cpx)
if COMPDIM == 3
for it = 1:numel(FRAMES)
tmp = squeeze(compr(OPE_.*FLD_(:,:,:,it)));
FIELD(:,:,it) = squeeze(process(tmp));
end
else
if REALP
tmp = zeros(Ny,Nx,Nz);
else
tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
end
for it = 1:numel(FRAMES)
for iz = 1:numel(DATA.z)
tmp(:,:,iz) = squeeze(process(OPE_.*FLD_(:,:,:,it)));
end
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
TOPLOT.FIELD = FIELD;
TOPLOT.FRAMES = FRAMES;
TOPLOT.X = shift_x(X);
TOPLOT.Y = shift_y(Y);
TOPLOT.FIELDNAME = FIELDNAME;
TOPLOT.XNAME = XNAME;
TOPLOT.YNAME = YNAME;
TOPLOT.FILENAME = [NAME,'_',OPTIONS.PLAN,'_',COMPNAME,'_',POLARNAME];
TOPLOT.DIMENSIONS= DIMENSIONS;
TOPLOT.ASPECT = ASPECT;
TOPLOT.FRAMES = FRAMES;
TOPLOT.INTERP = INTERP;
end
wk/Ajay_scan_CH4_lin_ITG.m
+ 122
− 66
View file @ 44d4d8b2
gyacomodir = '/home/ahoffman/gyacomo/';
gyacomodir = pwd;
gyacomodir = gyacomodir(1:end-2);
addpath(genpath([gyacomodir,'matlab'])) % ... add
addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
default_plots_options
% EXECNAME = 'gyacomo_dbg';
EXECNAME = 'gyacomo';
EXECNAME = 'gyacomo_alpha';
% EXECNAME = 'gyacomo';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
NU_a = [0:0.01:0.05]*0.1/0.045;
NU_a = [0 0.002 0.005 0.01 0.02 0.05]/0.45;
% P_a = [2 4 6 8 10 12 16];
% NU_a = 0.1;
P_a = [2 4];
% NU_a = 0.0;
% P_a = [4 8 12 20];
P_a = [20];
J_a = floor(P_a/2);
SIMID = 'Ajay_CH4_lin_ITG'; % Name of the simulation
RUN = 1; % to make sure it does not try to run
RERUN = 1; % rerun if the data does not exist
CO = 'DG';
GKCO = 0; % gyrokinetic operator
COLL_KCUT = 1.75;
% model
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-3; % electron plasma beta
K_Ti = 8.0;
K_Ni = 2.0;
K_Te = 8.0;
K_Ne = 2.0;
GEOMETRY= 'miller';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
% DT = 2e-3;
DT = 1e-3;
TMAX = 20;
TMAX = 25;
kymin = 0.35;
Ny = 2;
SIMID = 'Ajay_CH4_lin_ITG'; % Name of the simulation
RUN = 1;
g_ky = zeros(numel(NU_a),numel(P_a),Ny/2);
NY = 2;
g_ky = zeros(numel(NU_a),numel(P_a),NY/2);
g_avg= g_ky*0;
g_std= g_ky*0;
......@@ -33,67 +49,87 @@ for P = P_a
i = 1;
for NU = NU_a
%Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
TAU = 1.0; % e/i temperature ratio
K_Ni = 2.0; K_Ne = K_Ni;
K_Te = K_Ti; % Temperature '''
%% PHYSICAL PARAMETERS
TAU = 1.0; % e/i temperature ratio
% SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-3; % electron plasma beta
J = P/2;
% J = 2;
PMAXE = P; JMAXE = J;
PMAXI = P; JMAXI = J;
NX = 16; % real space x-gridpoints
NY = Ny; % '' y-gridpoints
LX = 142.9; % Size of the squared frequency domain
LY = 2*pi/kymin;
%% GRID PARAMETERS
% P = 20;
J = floor(P/2);
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
NX = 8; % real space x-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 2*pi/kymin; % Size of the squared frequency domain
NZ = 24; % number of perpendicular planes (parallel grid)
NPOL = 1; SG = 0;
% GEOMETRY= 's-alpha';
GEOMETRY= 'miller';
NPOL = 1;
SG = 0; % Staggered z grids option
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
%% GEOMETRY
% GEOMETRY= 's-alpha';
EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
SPS0D = 1; SPS2D = -1; SPS3D = 1;SPS5D= 1/5; SPSCP = 0;
PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = -1; % Sampling per time unit for 2D arrays
SPS3D = 1; % Sampling per time unit for 2D arrays
SPS5D = 1/2; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
GKCO = 1; % gyrokinetic operator
% Collision operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom
INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
W_PHI = 1; W_NA00 = 1;
W_DENS = 1; W_TEMP = 1;
W_NAPJ = 1; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
HD_CO = 0.0; % Hyper diffusivity cutoff ratio
MU = 0.0; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; N_HD = 4;
MU_Z = 2.0; MU_P = 0.0; %
MU_J = 0.0; LAMBDAD = 0.0;
MU_Y = MU; %
N_HD = 4;
MU_Z = 0.2; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
GRADB = 1.0;CURVB = 1.0;
%%-------------------------------------------------------------------------
% RUN
GRADB = 1.0;
CURVB = 1.0;
%% RUN
setup
if RUN
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 1 4 0; cd ../../../wk'])
% naming
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
% check if data exist to run if no data
data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
if ((RERUN || isempty(data.NU_EVOL)) && RUN)
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
end
% Load results
% Load results after trying to run
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
......@@ -102,17 +138,29 @@ for NU = NU_a
options.TRANGE = [0.5 1]*data.Ts3D(end);
options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
lg = compute_fluxtube_growth_rate(data,options);
[gmax, kmax] = max(lg.g_ky(:,end));
[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
% lg = compute_fluxtube_growth_rate(data,options);
% [gmax, kmax] = max(lg.g_ky(:,end));
% [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
% g_ky(i,j,:) = g_ky;
%
% g_avg(i,j,:) = lg.avg_g;
% g_std(i,j,:) = lg.std_g;
g_ky(i,j,:) = lg.avg_g;
[~,it1] = min(abs(data.Ts3D-0.5*data.Ts3D(end))); % start of the measurement time window
[~,it2] = min(abs(data.Ts3D-1.0*data.Ts3D(end))); % end of ...
field = 0;
field_t = 0;
for ik = 1:NY/2+1
field = squeeze(sum(abs(data.PHI),3)); % take the sum over z
field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
to_measure = log(field_t);
gr = polyfit(data.Ts3D(it1:it2),to_measure(it1:it2),1);
g_ky(i,j,ik) = gr(1);
end
[gmax, ikmax] = max(g_ky(i,j,:));
g_avg(i,j,:) = lg.avg_g;
g_std(i,j,:) = lg.std_g;
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data.ky(ikmax)); disp(msg);
i = i + 1;
end
......@@ -123,23 +171,27 @@ if 1
%% Study of the peak growth rate
figure
y_ = g_avg;
e_ = g_std;
y_ = g_ky;
e_ = 0.05;
% filter to noisy data
y_ = y_.*(y_-e_>0);
e_ = e_ .* (y_>0);
[y_,idx_] = max(g_avg,[],3);
for i = 1:numel(idx_)
e_ = g_std(:,:,idx_(i));
end
[y_,idx_a] = max(g_ky,[],3);
% for i = 1:numel(idx_)
% e_ = g_std(:,:,idx_(i));
% end
colors_ = summer(numel(NU_a));
colors_ = lines(numel(NU_a));
subplot(121)
for i = 1:numel(NU_a)
errorbar(P_a,y_(i,:),e_(i,:),...
'LineWidth',1.2,...
% errorbar(P_a,y_(i,:),e_(i,:),...
% 'LineWidth',1.2,...
% 'DisplayName',['$\nu=$',num2str(NU_a(i))],...
% 'color',colors_(i,:));
plot(P_a,y_(i,:),'s-',...
'LineWidth',2.0,...
'DisplayName',['$\nu=$',num2str(NU_a(i))],...
'color',colors_(i,:));
hold on;
......@@ -150,9 +202,13 @@ legend('show'); xlabel('$P$, $J=P/2$'); ylabel('$\gamma$');
colors_ = jet(numel(P_a));
subplot(122)
for j = 1:numel(P_a)
errorbar(NU_a,y_(:,j),e_(:,j),...
'LineWidth',1.2,...
'DisplayName',['(',num2str(P_a(j)),',',num2str(P_a(j)/2),')'],...
% errorbar(NU_a,y_(:,j),e_(:,j),...
% 'LineWidth',1.2,...
% 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
% 'color',colors_(j,:));
plot(NU_a,y_(:,j),'s-',...
'LineWidth',2.0,...
'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
'color',colors_(j,:));
hold on;
end
......
wk/CBC_nu_CO_scan_miller.m
+ 4
− 4
View file @ 44d4d8b2
......@@ -15,10 +15,10 @@ SIMID = 'convergence_pITG_miller'; % Name of the simulation
RUN = 0; % to make sure it does not try to run
RERUN = 0; % rerun if the data does not exist
NU_a = [0.02];
P_a = [24];
% NU_a = [0.0 0.01 0.02 0.05 0.1];
% P_a = [2 4:4:28];
% NU_a = [0.02];
% P_a = [24];
NU_a = [0.0 0.01 0.02 0.05 0.1];
P_a = [2 4:4:28];
% P_a = [24:4:28];
J_a = floor(P_a/2);
% collision setting
......
wk/analysis_gyacomo.m
+ 15
− 13
View file @ 44d4d8b2
......@@ -62,20 +62,20 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
options.NAME = '\phi';
% options.NAME = '\phi';
% options.NAME = '\omega_z';
% options.NAME = 'N_i^{00}';
% options.NAME = 'v_x';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
options.PLAN = 'yz';
options.NAME = 'n_i-n_e';
options.PLAN = 'xy';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
% options.TIME = data.Ts3D;
options.TIME = [1:1:9000];
options.TIME = [500:1:9000];
data.EPS = 0.1;
data.a = data.EPS * 2000;
options.RESOLUTION = 256;
......@@ -83,23 +83,25 @@ create_film(data,options,'.gif')
end
if 0
%% 2D snapshots
%% fields snapshots
% Options
options.INTERP = 0;
options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
options.NORMALIZE = 0;
options.NAME = '\phi';
% options.NAME = '\psi';
% options.NAME = '\omega_z';
% options.NAME = 'n_i';
% options.NAME = 'n_i-n_e';
% options.NAME = '\phi^{NZ}';
% options.NAME = 'N_i^{00}';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 's_{Ex}';
% options.NAME = 'Q_x';
% options.NAME = 'k^2n_e';
options.PLAN = 'kxz';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.PLAN = 'xy';
options.COMP = 'avg';
options.TIME = [100 250 300];
options.TIME = [800];
options.RESOLUTION = 256;
data.a = data.EPS * 2e3;
......@@ -108,7 +110,7 @@ fig = photomaton(data,options);
end
if 0
%% 3D plot on the geometry
%% plot on the geometry
options.INTERP = 0;
options.NAME = '\phi';
options.PLANES = [1];
......
wk/gene_data.FIGDIR deleted 100644 → 0
+ 0
− 0
View file @ 2bf9db7d
wk/header_2DZP_results.m
+ 2
− 2
View file @ 44d4d8b2
......@@ -180,7 +180,7 @@ resdir ='';
% resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6';
% resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6_conv';
% resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6_mu_0.5';
% resdir ='Zpinch_rerun/Kn_1.6_256x128x21x11';
resdir ='Zpinch_rerun/Kn_1.6_256x128x21x11';
%% nu scan
% resdir = 'Zpinch_rerun/kN_2.2_coll_scan_128x48x5x3';
......@@ -237,7 +237,7 @@ resdir ='';
% resdir = 'Zpinch_rerun/DGGK_kN_1.7_200x64x9x5_nu_0.01';
% resdir = 'Zpinch_rerun/LRGK_kN_2.4_300x98x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/LDGK_kN_1.9_200x64x5x3_nu_0.01';
resdir = 'Zpinch_rerun/COLLESS_kN_1.7_200x64x5x3';
% resdir = 'Zpinch_rerun/COLLESS_kN_1.7_200x64x5x3';
% resdir = 'Zpinch_rerun/COLLESS_kN_1.7_200x64x9x5';
JOBNUMMIN = 00; JOBNUMMAX = 10;
......
wk/lin_ITG.m
+ 11
− 11
View file @ 44d4d8b2
......@@ -11,7 +11,7 @@ addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
SIMID = 'dbg'; % Name of the simulation
RUN = 1; % To run or just to load
RUN = 0; % To run or just to load
default_plots_options
% EXECNAME = 'gyacomo_dbg';
EXECNAME = 'gyacomo_alpha';
......@@ -32,36 +32,36 @@ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-4; % electron plasma beta
%% GRID PARAMETERS
P = 8;
P = 20;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
NX = 8; % real space x-gridpoints
NX = 8; % real space x-gridpoints
NY = 2; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 2*pi/0.4; % Size of the squared frequency domain
LX = 2*pi/0.1; % Size of the squared frequency domain
LY = 2*pi/0.3; % Size of the squared frequency domain
NZ = 24; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
%% GEOMETRY
GEOMETRY= 's-alpha';
% GEOMETRY= 'miller';
% GEOMETRY= 's-alpha';
GEOMETRY= 'miller';
EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
% PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
TMAX = 50; % Maximal time unit
% DT = 1e-2; % Time step
DT = 1e-3; % Time step
DT = 5e-4; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = -1; % Sampling per time unit for 2D arrays
SPS3D = 1; % Sampling per time unit for 2D arrays
......@@ -112,7 +112,7 @@ setup
% Run linear simulation
if RUN
% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 2 2 1 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment