From 44d4d8b25bb23274820dfc47691ed2f7990f475d Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Tue, 10 Jan 2023 13:41:47 +0100
Subject: [PATCH] script save
---
matlab/load/quick_gene_load.m | 4 +-
matlab/plot/create_film.m | 2 +-
matlab/plot/photomaton.m | 11 +-
matlab/process_field_v2.m | 258 ++++++++++++++++++++++++++++++++++
wk/Ajay_scan_CH4_lin_ITG.m | 188 ++++++++++++++++---------
wk/CBC_nu_CO_scan_miller.m | 8 +-
wk/analysis_gyacomo.m | 28 ++--
wk/gene_data.FIGDIR | 0
wk/header_2DZP_results.m | 4 +-
wk/lin_ITG.m | 22 +--
10 files changed, 423 insertions(+), 102 deletions(-)
create mode 100644 matlab/process_field_v2.m
delete mode 100644 wk/gene_data.FIGDIR
diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index 55699d4e..4d1fc2bd 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -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';
diff --git a/matlab/plot/create_film.m b/matlab/plot/create_film.m
index 58b1baef..5f7aa15a 100644
--- a/matlab/plot/create_film.m
+++ b/matlab/plot/create_film.m
@@ -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];
diff --git a/matlab/plot/photomaton.m b/matlab/plot/photomaton.m
index c55123e0..e60fd8cc 100644
--- a/matlab/plot/photomaton.m
+++ b/matlab/plot/photomaton.m
@@ -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;
diff --git a/matlab/process_field_v2.m b/matlab/process_field_v2.m
new file mode 100644
index 00000000..4525a46a
--- /dev/null
+++ b/matlab/process_field_v2.m
@@ -0,0 +1,258 @@
+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
+
diff --git a/wk/Ajay_scan_CH4_lin_ITG.m b/wk/Ajay_scan_CH4_lin_ITG.m
index cb633fcf..1c1d2e5a 100644
--- a/wk/Ajay_scan_CH4_lin_ITG.m
+++ b/wk/Ajay_scan_CH4_lin_ITG.m
@@ -1,30 +1,46 @@
-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
diff --git a/wk/CBC_nu_CO_scan_miller.m b/wk/CBC_nu_CO_scan_miller.m
index 87b104ed..fcece647 100644
--- a/wk/CBC_nu_CO_scan_miller.m
+++ b/wk/CBC_nu_CO_scan_miller.m
@@ -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
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index 50b84ffa..423c3a09 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -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];
diff --git a/wk/gene_data.FIGDIR b/wk/gene_data.FIGDIR
deleted file mode 100644
index e69de29b..00000000
diff --git a/wk/header_2DZP_results.m b/wk/header_2DZP_results.m
index 6098397a..0b25b98b 100644
--- a/wk/header_2DZP_results.m
+++ b/wk/header_2DZP_results.m
@@ -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;
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index a54a321a..39cc2067 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -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'])
--
GitLab