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

Scripts update

parent b566ebf2
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 349 additions and 102 deletions
matlab/compile_results.m
+ 1
− 1
View file @ 53ffe66b
......@@ -277,7 +277,7 @@ else
% grids
DATA.Pe = Pe; DATA.Pi = Pi;
DATA.Je = Je; DATA.Ji = Ji;
DATA.kx = kx; DATA.ky = ky; DATA.z = z;
DATA.kx = kx; DATA.ky = ky; DATA.z = z; DATA.Npol = -z(1)/pi;
DATA.x = x; DATA.y = y;
DATA.ikx0 = ikx0; DATA.iky0 = iky0;
DATA.Nx = Nx; DATA.Ny = Ny; DATA.Nz = Nz; DATA.Nkx = Nkx; DATA.Nky = Nky;
......
matlab/compute/compute_3D_zpinch_growth_rate.m
+ 39
− 46
View file @ 53ffe66b
......@@ -4,7 +4,7 @@ FIGURE.FIGNAME = ['growth_rate_kx=0_ky=0_planes',DATA.PARAMS];
t = DATA.Ts3D;
[~,its] = min(abs(t-TRANGE(1)));
[~,ite] = min(abs(t-TRANGE(end)));
nkx = DATA.Nkx; nky = DATA.Nky; nkz = DATA.Nz;
nkx = DATA.Nkx; nky = DATA.Nky; nkz = DATA.Nz/2;
% Remove AA part
if DATA.Nx > 1
ikxnz = abs(DATA.PHI(1,:,1,1)) > 0;
......@@ -15,84 +15,77 @@ ikynz = abs(DATA.PHI(:,1,1,1)) > 0;
ikynz(1) = 1; ikxnz(1) = 1; %put k=0 in the analysis
phi = fft(DATA.PHI(ikynz,ikxnz,:,:),[],3);
Y = fft(DATA.PHI(ikynz,ikxnz,:,:),[],3);
phi = Y(:,:,2:2:end,:);
omega = zeros(nky,nkx,nkz);
for iz = 1:nkz
for iy = 1:nky
for ikz = 1:nkz
for iky = 1:nky
for ix = 1:nkx
omega(iy,ix,iz) = LinearFit_s(t(its:ite),squeeze(abs(phi(iy,ix,iz,its:ite))));
% omega(iy,ix,iz) = LinearFit_s(t(its:ite),squeeze(abs(phi(iy,ix,iz,its:ite))));
to_measure = squeeze(log(phi(iky,ix,ikz,its:ite)));
tmp = polyfit(t(its:ite),to_measure(:),1);
if ~(isnan(tmp(1)) || isinf(tmp(1)))
omega(iky,ix,ikz) = tmp(1);
end
end
end
end
%% plot
kx = DATA.kx; ky = DATA.ky; kz = [(0:nkz/2), (-nkz/2+1):-1];
poskx = kx>=0;
posky = ky>=0;
poskz = kz>=0;
kx = DATA.kx; ky = DATA.ky; kz = [(0:nkz/2), (-nkz/2+1):-1]/DATA.Npol;
kxeq0 = kx==0;
kyeq0 = ky==0;
kzeq0 = kz==0;
omega = omega(posky,poskx,poskz);
% kx = fftshift(kx,1);
% ky = fftshift(ky,1);
% kz = fftshift(kz,1);
FIGURE.fig = figure;
nplots = OPTIONS.kxky + OPTIONS.kzkx + OPTIONS.kzky;
iplot = 1;
if OPTIONS.kxky
[Y_XY,X_XY] = meshgrid(ky(posky),kx(poskx));
[Y_XY,X_XY] = meshgrid(ky,kx);
subplot(1,nplots,iplot)
if ~OPTIONS.keq0
toplot = squeeze(max(real(omega(:,:,:)),[],3));
pclr= pcolor(X_XY,Y_XY,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_x$'); ylabel('$k_y$');
title('$\max(\gamma)_{kz}$');
else
toplot = squeeze(real(omega(:,:,kzeq0)));
pclr= pcolor(X_XY,Y_XY,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_x$'); ylabel('$k_y$');
title('$\gamma(k_z=0)$');
end
toplot = squeeze(real(omega(:,:,kzeq0)));
toplot = fftshift(toplot,2);
X_XY = fftshift(X_XY,1);
Y_XY = fftshift(Y_XY,1);
pclr= pcolor(X_XY,Y_XY,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_x$'); ylabel('$k_y$');
title('$\gamma(k_z=0)$');
if OPTIONS.INTERP; shading interp; end
caxis(max(max(abs(toplot))).*[-1,1]);
iplot = iplot + 1;
end
if OPTIONS.kzky
[Y_ZY,Z_ZY] = meshgrid(ky(posky),kz(poskz));
[Y_ZY,Z_ZY] = meshgrid(ky,kz);
subplot(1,nplots,iplot)
if ~OPTIONS.keq0
toplot = squeeze(max(real(omega(:,:,:)),[],1));
pclr= pcolor(Z_ZY,Y_ZY,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_x$'); ylabel('$k_y$');
title('$\max(\gamma)_{kx}$');
else
toplot = squeeze(real(omega(:,kxeq0,:)));
pclr= pcolor(Z_ZY,Y_ZY,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_z$'); ylabel('$k_y$');
title('$\gamma(k_x=0)$');
end
toplot = squeeze(real(omega(:,kxeq0,:)));
toplot = fftshift(toplot,2);
Z_ZY = fftshift(Z_ZY,1);
Y_ZY = fftshift(Y_ZY,1);
pclr= pcolor(Z_ZY,Y_ZY,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_z$'); ylabel('$k_y$');
title('$\gamma(k_x=0)$');
if OPTIONS.INTERP; shading interp; end
caxis(max(max(abs(toplot))).*[-1,1]);
iplot = iplot + 1;
end
if OPTIONS.kzkx
[X_ZX,Z_ZX] = meshgrid(kx(poskx),kz(poskz));
[X_ZX,Z_ZX] = meshgrid(kx,kz);
subplot(1,nplots,iplot)
if ~OPTIONS.keq0
toplot = squeeze(max(real(omega(:,:,:)),[],2));
pclr= pcolor(Z_ZX,X_ZX,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_z$'); ylabel('$k_x$');
title('$\max(\gamma)_{ky}$');
else
toplot = squeeze(real(omega(kyeq0,:,:)));
pclr= pcolor(Z_ZX,X_ZX,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_z$'); ylabel('$k_x$');
title('$\gamma(k_y=0)$');
end
toplot = squeeze(real(omega(kyeq0,:,:)));
toplot = fftshift(toplot,2);
Z_ZX = fftshift(Z_ZX,1);
X_ZX = fftshift(X_ZX,1);
pclr= pcolor(Z_ZX,X_ZX,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_z$'); ylabel('$k_x$');
title('$\gamma(k_y=0)$');
end
if OPTIONS.INTERP; shading interp; end
caxis(max(max(abs(toplot))).*[-1,1]);
......
matlab/load/load_params.m
+ 6
− 0
View file @ 53ffe66b
......@@ -24,6 +24,12 @@ catch
DATA.K_T = h5readatt(filename,'/data/input','k_Ti');
end
end
DATA.sigma_e = h5readatt(filename,'/data/input','sigma_e');
DATA.sigma_i = h5readatt(filename,'/data/input','sigma_i');
DATA.tau_e = h5readatt(filename,'/data/input','tau_e');
DATA.tau_i = h5readatt(filename,'/data/input','tau_i');
DATA.q_e = h5readatt(filename,'/data/input','q_e');
DATA.q_i = h5readatt(filename,'/data/input','q_i');
DATA.Q0 = h5readatt(filename,'/data/input','q0');
DATA.SHEAR = h5readatt(filename,'/data/input','shear');
DATA.EPS = h5readatt(filename,'/data/input','eps');
......
matlab/plot/show_moments_spectrum.m
+ 4
− 4
View file @ 53ffe66b
......@@ -22,7 +22,7 @@ if ~OPTIONS.ST
[~,it0] = min(abs(t0-DATA.Ts5D));
[~,it1] = min(abs(t1-DATA.Ts5D));
Nipj = mean(Nipj(:,:,it0:it1),3);
if DATA.K_E
if DATA.KIN_E
Nepj = mean(Nepj(:,:,it0:it1),3);
end
if numel(OPTIONS.TIME) == 1
......@@ -34,13 +34,13 @@ if ~OPTIONS.ST
ymini = min(Nipj); ymaxi = max(Nipj);
if DATA.K_E
if DATA.KIN_E
Nepj = squeeze(Nepj);
ymine = min(Nepj); ymaxe = max(Nepj);
ymax = max([ymaxi ymaxe]);
ymin = min([ymini ymine]);
end
if DATA.K_E
if DATA.KIN_E
subplot(121)
end
if ~OPTIONS.P2J
......@@ -60,7 +60,7 @@ if ~OPTIONS.ST
legend('show');
title([TITLE,' He-La ion spectrum']);
if DATA.K_E
if DATA.KIN_E
subplot(122)
if ~OPTIONS.P2J
for ij = 1:numel(Je)
......
matlab/process_field_v2.m
+ 31
− 3
View file @ 53ffe66b
......@@ -65,7 +65,29 @@ else
FIELD = zeros(sz(1),sz(2),Nt);
end
%% Process the field to plot
% short term writing --
b_e = DATA.sigma_e*sqrt(2*DATA.tau_e)*sqrt(KX.^2+KY.^2);
adiab_e = kernel(0,b_e);
pol_e = kernel(0,b_e).^2;
for n = 1:DATA.Jmaxe
adiab_e = adiab_e + kernel(n,b_e).^2;
pol_e = pol_e + kernel(n,b_e).^2;
end
adiab_e = DATA.q_e/DATA.tau_e .* adiab_e;
pol_e = DATA.q_e^2/DATA.tau_e * (1 - pol_e);
b_i = DATA.sigma_i*sqrt(2*DATA.tau_i)*sqrt(KX.^2+KY.^2);
adiab_i = kernel(0,b_i);
pol_i = kernel(0,b_i).^2;
for n = 1:DATA.Jmaxi
adiab_i = adiab_i + kernel(n,b_i).^2;
pol_i = pol_i + kernel(n,b_i).^2;
end
pol_i = DATA.q_i^2/DATA.tau_i * (1 - pol_i);
adiab_i = DATA.q_i/DATA.tau_i .* adiab_i;
poisson_op = (pol_i + pol_e);
% --
switch OPTIONS.COMP
case 'avg'
compr = @(x) mean(x,COMPDIM);
......@@ -130,6 +152,7 @@ end
switch OPTIONS.NAME
case '\phi' %ES pot
NAME = 'phi';
% FLD_ = DATA.PHI(:,:,:,FRAMES); % data to plot
FLD_ = DATA.PHI(:,:,:,FRAMES); % data to plot
OPE_ = 1; % Operation on data
case '\psi' %EM pot
......@@ -168,13 +191,17 @@ switch OPTIONS.NAME
NAME = 'Ne00';
FLD_ = DATA.Ne00(:,:,:,FRAMES);
OPE_ = 1;
case 'N_i^{00}-N_e^{00}' % first electron gyromoment
NAME = 'Ni00-Ne00';
FLD_ = (DATA.Ni00(:,:,:,FRAMES)-DATA.Ne00(:,:,:,FRAMES))./(poisson_op+(poisson_op==0));
OPE_ = 1;
case 'n_i' % ion density
NAME = 'ni';
FLD_ = DATA.DENS_I(:,:,:,FRAMES);
FLD_ = DATA.DENS_I(:,:,:,FRAMES) - adiab_i.* DATA.PHI(:,:,:,FRAMES);
OPE_ = 1;
case 'n_e' % electron density
NAME = 'ne';
FLD_ = DATA.DENS_E(:,:,:,FRAMES);
FLD_ = DATA.DENS_E(:,:,:,FRAMES) - adiab_e.* DATA.PHI(:,:,:,FRAMES);
OPE_ = 1;
case 'k^2n_e' % electron vorticity
NAME = 'k2ne';
......@@ -183,7 +210,8 @@ switch OPTIONS.NAME
case 'n_i-n_e' % polarisation
NAME = 'pol';
OPE_ = 1;
FLD_ = DATA.DENS_I(:,:,:,FRAMES)+DATA.DENS_E(:,:,:,FRAMES);
FLD_ = ((DATA.DENS_I(:,:,:,FRAMES)- adiab_i.* DATA.PHI(:,:,:,FRAMES))...
-(DATA.DENS_E(:,:,:,FRAMES)- adiab_e.* DATA.PHI(:,:,:,FRAMES)));
case 'T_i' % ion temperature
NAME = 'Ti';
FLD_ = DATA.TEMP_I(:,:,:,FRAMES);
......
wk/analysis_gyacomo.m
+ 7
− 5
View file @ 53ffe66b
......@@ -14,7 +14,7 @@ MISCDIR = ['/misc/gyacomo_outputs/',resdir,'/']; %For long term storage
system(['mkdir -p ',MISCDIR]);
system(['mkdir -p ',LOCALDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
system(CMD);
% system(CMD);
% Load outputs from jobnummin up to jobnummax
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
data.localdir = LOCALDIR;
......@@ -68,14 +68,15 @@ options.POLARPLOT = 0;
% options.NAME = 'v_x';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
options.NAME = 'n_i-n_e';
options.NAME = 'n_e';
% 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 = [500:1:9000];
options.TIME = [00:1:9000];
data.EPS = 0.1;
data.a = data.EPS * 2000;
options.RESOLUTION = 256;
......@@ -92,16 +93,17 @@ options.NORMALIZE = 0;
options.NAME = '\phi';
% options.NAME = '\psi';
% options.NAME = '\omega_z';
% options.NAME = 'n_i';
% options.NAME = 'n_e';
% options.NAME = 'n_i-n_e';
% options.NAME = '\phi^{NZ}';
% options.NAME = 'N_i^{00}';
% options.NAME = 'N_i^{00}-N_e^{00}';
% options.NAME = 's_{Ex}';
% options.NAME = 'Q_x';
% options.NAME = 'k^2n_e';
options.PLAN = 'xy';
options.COMP = 'avg';
options.TIME = [800];
options.TIME = [500];
options.RESOLUTION = 256;
data.a = data.EPS * 2e3;
......
wk/header_2DZP_results.m
+ 3
− 3
View file @ 53ffe66b
......@@ -194,7 +194,7 @@ resdir ='Zpinch_rerun/Kn_1.6_256x128x21x11';
% resdir ='Zpinch_rerun/nu_0.1_FCGK_200x48x5x3_kN_scan';
% resdir ='Zpinch_rerun/nu_0.1_SGGK_200x48x5x3_kN_scan';
% resdir = 'Zpinch_rerun/kN_1.7_FCGK_200x32x5x3_nu_scan';
% % resdir = 'Zpinch_rerun/kN_1.7_SGGK_200x32x7x4_nu_scan';
% resdir = 'Zpinch_rerun/kN_1.7_SGGK_200x32x7x4_nu_scan';
% resdir = 'Zpinch_rerun/kN_2.2_SGGK_200x32x5x3_nu_scan';
% resdir = 'Zpinch_rerun/kN_1.7_SGGK_256x64x5x3_nu_scan';
......@@ -230,7 +230,7 @@ resdir ='Zpinch_rerun/Kn_1.6_256x128x21x11';
% resdir = 'Zpinch_rerun/DGGK_kN_scan_200x64x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/even_DGGK_kN_scan_200x64x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/LRGK_kN_scan_200x64x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/LDGK_kN_scan_200x64x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/LDGK_kN_scan_202x64x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/DGGK_kN_1.6_200x64x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/DGGK_kN_1.7_200x64x5x3_nu_0.01';
......@@ -240,6 +240,6 @@ resdir ='Zpinch_rerun/Kn_1.6_256x128x21x11';
% resdir = 'Zpinch_rerun/COLLESS_kN_1.7_200x64x5x3';
% resdir = 'Zpinch_rerun/COLLESS_kN_1.7_200x64x9x5';
JOBNUMMIN = 00; JOBNUMMAX = 10;
JOBNUMMIN = 01; JOBNUMMAX = 03;
resdir = ['results/',resdir];
run analysis_gyacomo
\ No newline at end of file
wk/lin_3Dzpinch.m 0 → 100644
+ 201
− 0
View file @ 53ffe66b
%% QUICK RUN SCRIPT for linear entropy mode (ETPY) in a Zpinch
% This script create a directory in /results and run a simulation directly
% from matlab framework. It is meant to run only small problems in linear
% for benchmark and debugging purpose since it makes matlab "busy"
%
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';
SIMID = 'dbg'; % Name of the simulation
RUN = 0; % To run or just to load
default_plots_options
% EXECNAME = 'gyacomo_debug';
% EXECNAME = 'gyacomo_alpha';
EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
TAU = 1e-2; % e/i temperature ratio
K_Ne = 0; % ele Density '''
K_Te = 0; % ele Temperature '''
K_Ni = 0; % ion Density gradient drive
K_Ti = 70; % ion Temperature '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 0.0; % electron plasma beta
%% GRID PARAMETERS
P = 4;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
NX = 20; % real space x-gridpoints
NY = 20; % '' y-gridpoints
LX = 2*pi/0.4; % Size of the squared frequency domain
LY = 2*pi/0.4; % Size of the squared frequency domain
NZ = 48; % 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
GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
SHEAR = 0.0; % 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'
SHIFT_Y = 0.0;
%% TIME PARMETERS
TMAX = 50; % Maximal time unit
DT = 5e-2; % 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
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)
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
CO = 'DG';
GKCO = 0; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
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 = 1.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
GRADB = 0.1;
CURVB = 0.1;
COLL_KCUT = 1000;
%%-------------------------------------------------------------------------
%% RUN
setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
% 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 6 ',gyacomodir,'bin/',EXECNAME,' 1 3 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
end
%% Load results
%%
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
FIGDIR = LOCALDIR;
% Load outputs from jobnummin up to jobnummax
JOBNUMMIN = 00; JOBNUMMAX = 01;
data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
data.FIGDIR = ['../results/',SIMID,'/',PARAMS,'/'];
%% Short analysis
if 0
%% linear growth rate (adapted for 2D zpinch and fluxtube)
options.TRANGE = [0.5 1]*data.Ts3D(end);
options.NPLOTS = 3; % 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);
end
if 0
%% Ballooning plot
options.time_2_plot = [10 50];
options.kymodes = [0.1 0.2 0.4];
options.normalized = 1;
% options.field = 'phi';
fig = plot_ballooning(data,options);
end
if 0
%% Hermite-Laguerre spectrum
% options.TIME = 'avg';
options.P2J = 0;
options.ST = 0;
options.PLOT_TYPE = 'space-time';
% options.PLOT_TYPE = 'Tavg-1D';ls
% options.PLOT_TYPE = 'Tavg-2D';
options.NORMALIZED = 1;
options.JOBNUM = 0;
options.TIME = [0 50];
options.specie = 'i';
options.compz = 'avg';
fig = show_moments_spectrum(data,options);
% fig = show_napjz(data,options);
% save_figure(data,fig)
end
if 1
%% linear growth rate for 3D Zpinch (kz fourier transform)
trange = [0.5 1]*data.Ts3D(end);
options.INTERP = 0;
options.kxky = 0;
options.kzkx = 0;
options.kzky = 1;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
end
if 0
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 1;
options.TIME = [0:1000];
% options.TIME = [0.5:0.01:1].*data.Ts3D(end);
options.NMA = 1;
options.NMODES = 32;
options.iz = 'avg';
options.ik = 1;
fig = mode_growth_meter(data,options);
save_figure(data,fig,'.png')
end
if 0
%% RH TEST
ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end);
[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));
plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3));
figure
plot(data.Ts3D(it0:it1), plt(data.PHI));
xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky)))
end
\ No newline at end of file
wk/lin_ETPY.m
+ 51
− 34
View file @ 53ffe66b
......@@ -3,13 +3,17 @@
% from matlab framework. It is meant to run only small problems in linear
% for benchmark and debugging purpose since it makes matlab "busy"
%
SIMID = 'lin_ETPY'; % Name of the simulation
% SIMID = 'dbg'; % Name of the simulation
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';
SIMID = 'dbg'; % Name of the simulation
RUN = 1; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
PROGDIR = '/home/ahoffman/gyacomo/';
EXECNAME = 'gyacomo_dbg';
% EXECNAME = 'gyacomo_debug';
EXECNAME = 'gyacomo_alpha';
% EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
......@@ -26,35 +30,40 @@ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 0.0; % electron plasma beta
%% GRID PARAMETERS
P = 8;
P = 4;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
NX = 2; % real space x-gridpoints
NY = 100; % '' y-gridpoints
NY = 40; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 200;%2*pi/0.05; % Size of the squared frequency domain
LY = 2*pi/0.05; % Size of the squared frequency domain
NZ = 1; % 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
GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
KAPPA = 0.0; % elongation
SHEAR = 0.0; % 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))
EPS = 0.18; % inverse aspect ratio
PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
TMAX = 500; % Maximal time unit
TMAX = 50; % Maximal time unit
% DT = 1e-2; % Time step
DT = 1e-2; % 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
SPS5D = 1/5; % Sampling per time unit for 5D arrays
SPS5D = 1/2; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
......@@ -62,14 +71,14 @@ LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
CO = 'DG';
GKCO = 1; % gyrokinetic operator
GKCO = 0; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom
NUMERICAL_SCHEME = 'RK2'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
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;
......@@ -85,7 +94,7 @@ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
MU_Z = 2.0; %
MU_Z = 0.2; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
......@@ -100,23 +109,27 @@ setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',PROGDIR,'bin/',EXECNAME,' 1 1 1 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,'/; time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
end
%% Load results
%%
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [PROGDIR,'results/',filename,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
FIGDIR = LOCALDIR;
% Load outputs from jobnummin up to jobnummax
JOBNUMMIN = 00; JOBNUMMAX = 00;
JOBNUMMIN = 00; JOBNUMMAX = 01;
data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
data.FIGDIR = ['../results/',SIMID,'/',PARAMS,'/'];
%% Short analysis
if 1
if 0
%% linear growth rate (adapted for 2D zpinch and fluxtube)
options.TRANGE = [0.5 1]*data.Ts3D(end);
options.NPLOTS = 2; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
options.NPLOTS = 3; % 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));
......@@ -127,8 +140,8 @@ end
if 0
%% Ballooning plot
options.time_2_plot = [120];
options.kymodes = [0.9];
options.time_2_plot = [10 50];
options.kymodes = [0.1 0.2 0.4];
options.normalized = 1;
% options.field = 'phi';
fig = plot_ballooning(data,options);
......@@ -140,9 +153,10 @@ if 0
options.P2J = 1;
options.ST = 1;
options.PLOT_TYPE = 'space-time';
% options.PLOT_TYPE = 'Tavg-1D';
% options.PLOT_TYPE = 'Tavg-1D';ls
% options.PLOT_TYPE = 'Tavg-2D';
options.NORMALIZED = 0;
options.NORMALIZED = 1;
options.JOBNUM = 0;
options.TIME = [0 50];
options.specie = 'i';
......@@ -155,21 +169,24 @@ end
if 0
%% linear growth rate for 3D Zpinch (kz fourier transform)
trange = [0.5 1]*data.Ts3D(end);
options.INTERP = 0;
options.keq0 = 1; % chose to plot planes at k=0 or max
options.kxky = 1;
options.kzkx = 0;
options.kzky = 0;
options.kzkx = 1;
options.kzky = 1;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
save_figure(data,fig)
end
if 0
if 1
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 10;
options.K2PLOT = 1;
options.TIME = [0:1000];
% options.TIME = [0.5:0.01:1].*data.Ts3D(end);
options.NMA = 1;
options.NMODES = 10;
options.NMODES = 32;
options.iz = 'avg';
options.ik = 1;
fig = mode_growth_meter(data,options);
save_figure(data,fig,'.png')
end
......@@ -184,4 +201,4 @@ figure
plot(data.Ts3D(it0:it1), plt(data.PHI));
xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky)))
end
end
\ No newline at end of file
wk/lin_ITG.m
+ 6
− 6
View file @ 53ffe66b
......@@ -11,9 +11,9 @@ 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 = 0; % To run or just to load
RUN = 1; % To run or just to load
default_plots_options
% EXECNAME = 'gyacomo_dbg';
% EXECNAME = 'gyacomo_debug';
EXECNAME = 'gyacomo_alpha';
% EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......@@ -29,10 +29,10 @@ K_Ni = 1*2.22; % ion Density gradient drive
K_Ti = 6.96; % ion Temperature '''
% SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-4; % electron plasma beta
%% GRID PARAMETERS
P = 20;
P = 4;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
......@@ -60,8 +60,8 @@ PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
TMAX = 50; % Maximal time unit
% DT = 1e-2; % Time step
DT = 5e-4; % Time step
DT = 1e-2; % 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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment