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

upgrade scripts

parent 69d071f6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
matlab/compute/compute_fluxtube_growth_rate.m
+ 25
− 16
View file @ c35f0204
function [ linear_gr ] = compute_fluxtube_growth_rate(DATA, TRANGE, PLOT)
function [ linear_gr ] = compute_fluxtube_growth_rate(DATA, OPTIONS)
% Remove AA part
if DATA.Nx > 1
......@@ -11,8 +11,8 @@ ikynz = logical(ikynz .* (DATA.ky > 0));
phi = DATA.PHI(ikynz,ikxnz,:,:);
t = DATA.Ts3D;
[~,its] = min(abs(t-TRANGE(1)));
[~,ite] = min(abs(t-TRANGE(end)));
[~,its] = min(abs(t-OPTIONS.TRANGE(1)));
[~,ite] = min(abs(t-OPTIONS.TRANGE(end)));
w_ky = zeros(sum(ikynz),ite-its);
ce = zeros(sum(ikynz),ite-its);
......@@ -34,7 +34,7 @@ for it = its+1:ite
end
[kys, Is] = sort(DATA.ky(ikynz));
linear_gr.trange = t(its:ite);
linear_gr.OPTIONS.TRANGE = t(its:ite);
linear_gr.g_ky = real(w_ky(Is,:));
linear_gr.w_ky = imag(w_ky(Is,:));
linear_gr.ce = abs(ce);
......@@ -44,20 +44,29 @@ linear_gr.avg_g = mean(real(w_ky(Is,:)),2);
linear_gr.std_w = std (imag(w_ky(Is,:)),[],2);
linear_gr.avg_w = mean(imag(w_ky(Is,:)),2);
if PLOT >0
if OPTIONS.NPLOTS >0
figure
if PLOT > 1
if OPTIONS.NPLOTS > 1
subplot(1,2,1)
end
x_ = linear_gr.ky;
plt = @(x) x;
OVERK = '';
if OPTIONS.GOK == 1
plt = @(x) x./x_;
OVERK = '/$k_y$';
elseif OPTIONS.GOK == 2
plt = @(x) x.^2./x_.^3;
OVERK = '/$k_y$';
end
plot(x_,plt(linear_gr.g_ky(:,end)),'-o','DisplayName',['$Re(\omega_{k_y})$',OVERK]); hold on;
plot(x_,plt(linear_gr.w_ky(:,end)),'--*','DisplayName',['$Im(\omega_{k_y})$',OVERK]); hold on;
plot(x_,plt(linear_gr.ce (:,end)),'x','DisplayName',['$\epsilon$',OVERK]); hold on;
plot(linear_gr.ky,linear_gr.g_ky(:,end),'-o','DisplayName','$Re(\omega_{k_y})$'); hold on;
plot(linear_gr.ky,linear_gr.w_ky(:,end),'--*','DisplayName','$Im(\omega_{k_y})$'); hold on;
plot(linear_gr.ky,linear_gr.ce (:,end),'x','DisplayName','$\epsilon$'); hold on;
errorbar(linear_gr.ky,linear_gr.avg_g,linear_gr.std_g,...
errorbar(linear_gr.ky,plt(linear_gr.avg_g),plt(linear_gr.std_g),...
'-o','DisplayName','$Re(\omega_{k_y})$',...
'LineWidth',1.5); hold on;
errorbar(linear_gr.ky,linear_gr.avg_w,linear_gr.std_w,...
errorbar(linear_gr.ky,plt(linear_gr.avg_w),plt(linear_gr.std_w),...
'--*','DisplayName','$Im(\omega_{k_y})$',...
'LineWidth',1.5); hold on;
......@@ -66,10 +75,10 @@ end
legend('show');
title(DATA.param_title);
if PLOT > 1
if PLOT == 2
if OPTIONS.NPLOTS > 1
if OPTIONS.NPLOTS == 2
subplot(1,2,2)
elseif PLOT == 3
elseif OPTIONS.NPLOTS == 3
subplot(2,2,2)
end
plot(DATA.Ts3D(its+1:ite),linear_gr.g_ky(Is,:)); hold on;
......@@ -77,7 +86,7 @@ if PLOT > 1
xlabel('t'); ylabel('$\gamma(t),\omega(t)$'); xlim([DATA.Ts3D(1) DATA.Ts3D(end)]);
end
if PLOT > 2
if OPTIONS.NPLOTS > 2
xlabel([]); xticks([]);
subplot(2,2,4)
[~,ikx0] = min(abs(DATA.kx));
......
matlab/evolve_tracers.m
+ 4
− 4
View file @ c35f0204
% Options
SHOW_FILM = 1;
SHOW_FILM = 0;
field2plot ='phi';
INIT = 'lin'; % lin (for a line)/ round (for a small round)/ gauss for random
U_TIME = 15000; % >0 for frozen velocity at a given time, -1 for evolving field
U_TIME = 1000; % >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 = 100;
dt_ = 0.1;
Nstep = ceil(Tfin/dt_);
% Init tracers
Np = 20; %number of tracers
Np = 200; %number of tracers
% color = tcolors;
color = jet(Np);
tcolors = distinguishable_colors(Np); %Their colors
......@@ -182,7 +182,7 @@ while(t_<Tfin && it <= Nstep)
% push the particle
% q = sign(-u___(3));
q = -u___(3);
q = 1;%-u___(3);
% q =1;
x_ = x_ + dt_*u___(1)*q;
y_ = y_ + dt_*u___(2)*q;
......
matlab/load/load_3D_data.m
+ 18
− 3
View file @ c35f0204
......@@ -8,9 +8,24 @@ function [ data, time, dt ] = load_3D_data( filename, variablename )
cstart= h5readatt(filename,'/data/input','start_iframe3d');
data = zeros(numel(ky),numel(kx),numel(z),numel(time));
for it = 1:numel(time)
tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
[KX,KY] = meshgrid(kx,ky);
switch variablename
case 'vEx'
for it = 1:numel(time)
tmp = h5read(filename,['/data/var3d/phi/', num2str(cstart+it,'%06d')]);
data(:,:,:,it) = 1i.*KY.*(tmp.real + 1i * tmp.imaginary);
end
case 'vEy'
for it = 1:numel(time)
tmp = h5read(filename,['/data/var3d/phi/', num2str(cstart+it,'%06d')]);
data(:,:,:,it) = -1i.*KX.*(tmp.real + 1i * tmp.imaginary);
end
otherwise
for it = 1:numel(time)
tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
end
end
end
......
wk/header_2DZP_results.m
+ 6
− 2
View file @ c35f0204
......@@ -195,10 +195,14 @@ resdir ='';
% 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_2.2_SGGK_200x32x5x3_nu_scan';
% resdir = 'Zpinch_rerun/kN_2.2_SGGK_200x32x5x3_nu_scan';
% resdir = 'Zpinch_rerun/kN_1.7_SGGK_256x64x5x3_nu_scan';
%% Convergence cases kN = (1.6 2.2) nu = (0.01 1.0)
resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x3x2_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x3x2_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x5x3_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x3x2_mu_0.01';
resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x7x4_mu_0.01';
JOBNUMMIN = 00; JOBNUMMAX = 10;
resdir = ['results/',resdir];
......
wk/lin_EPY.m
+ 14
− 10
View file @ c35f0204
......@@ -4,17 +4,19 @@
% for benchmark and debugging purpose since it makes matlab "busy"
%
SIMID = 'lin_EPY'; % Name of the simulation
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';
EXECNAME = 'gyacomo_dbg';
% EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
NU = 0.01; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_Ne = 2.2; % ele Density '''
K_Te = K_Ne/4; % ele Temperature '''
......@@ -24,14 +26,14 @@ 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 = 4;
P = 2;
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 = 2; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 120;%2*pi/0.05; % Size of the squared frequency domain
NZ = 1; % number of perpendicular planes (parallel grid)
......@@ -44,7 +46,7 @@ SHEAR = 0.8; % magnetic shear
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
TMAX = 50; % Maximal time unit
TMAX = 200; % Maximal time unit
DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
......@@ -56,7 +58,7 @@ JOB2LOAD= -1;
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
CO = 'LD';
CO = 'SG';
GKCO = 1; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
......@@ -93,7 +95,8 @@ 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 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'])
end
%% Load results
......@@ -107,9 +110,10 @@ data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from
%% Short analysis
if 1
%% linear growth rate (adapted for 2D zpinch and fluxtube)
trange = [0.5 1]*data.Ts3D(end);
nplots = 1; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
lg = compute_fluxtube_growth_rate(data,trange,nplots);
options.TRANGE = [0.5 1]*data.Ts3D(end);
options.NPLOTS = 1; % 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);
......
wk/lin_ITG.m
+ 1
− 1
View file @ c35f0204
......@@ -4,7 +4,7 @@
% for benchmark and debugging purpose since it makes matlab "busy"
%
SIMID = 'lin_ITG'; % Name of the simulation
RUN = 1; % To run or just to load
RUN = 0; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
HELAZDIR = '/home/ahoffman/gyacomo/';
......
wk/save_iFFT.m
+ 1
− 1
View file @ c35f0204
......@@ -3,7 +3,7 @@ simdir = '/misc/gyacomo_outputs/results/Zpinch_rerun';
resolu = 'UHD_512x256x2x1';
output = 'outputs_01.h5';
filename = [simdir,'/',resolu,'/',output];
fieldname = 'Ni00';
fieldname = 'vEy';
%% OUTPUT
outdir = '.';
......
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