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

update matlab

parent e829c731
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
matlab/profiler.m
+ 30
− 32
View file @ c36a0067
function [] = profiler(data)
function [] = profiler(DATADIR,JID)
%% load profiling
% filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00);
% outfilename = ['/misc/HeLaZ_outputs',filename(3:end)];
CPUTI=[]; DTSIM=[]; RHSTC=[]; POITC=[]; SAPTC=[]; COLTC=[];
GRATC=[]; NADTC=[]; ADVTC=[]; GHOTC=[]; CLOTC=[]; CHKTC=[];
DIATC=[]; EXBTC=[]; STETC=[]; TS0TC=[];
for i = 1:numel(data.outfilenames)
outfilename = data.outfilenames{i};
% for i = 1:numel(data.outfilenames)
outfilename = [DATADIR,'/',sprintf('outputs_%02d.h5',JID)];
CPUTI = [ CPUTI; double(h5readatt(outfilename,'/data/input','cpu_time'))];
DTSIM = [ DTSIM; h5readatt(outfilename,'/data/input/basic','dt')];
RHSTC = [ RHSTC; h5read(outfilename,'/profiler/Tc_rhs')];
POITC = [ POITC; h5read(outfilename,'/profiler/Tc_poisson')];
SAPTC = [ SAPTC; h5read(outfilename,'/profiler/Tc_Sapj')];
try
EXBTC = [ EXBTC; h5read(outfilename,'/profiler/Tc_ExBshear')];
catch
EXBTC = 0.*SAPTC;
end
COLTC = [ COLTC; h5read(outfilename,'/profiler/Tc_coll')];
GRATC = [ GRATC; h5read(outfilename,'/profiler/Tc_grad')];
NADTC = [ NADTC; h5read(outfilename,'/profiler/Tc_nadiab')];
......@@ -28,18 +24,18 @@ for i = 1:numel(data.outfilenames)
DIATC = [ DIATC; h5read(outfilename,'/profiler/Tc_diag')];
STETC = [ STETC; h5read(outfilename,'/profiler/Tc_step')];
TS0TC = [ TS0TC; h5read(outfilename,'/profiler/time')];
end
CPUTI = CPUTI(end);
% end
CPUTI = sum(CPUTI);
DTSIM = mean(DTSIM);
N_T = 13;
missing_Tc = STETC - RHSTC - ADVTC - GHOTC - CLOTC - EXBTC ...
MISTC = STETC - RHSTC - ADVTC - GHOTC - CLOTC - EXBTC ...
-COLTC - POITC - SAPTC - CHKTC - DIATC - GRATC - NADTC;
total_Tc = STETC;
TIME_PER_FCT = [diff(RHSTC); diff(ADVTC); diff(GHOTC);...
diff(CLOTC); diff(COLTC); diff(POITC); diff(SAPTC); diff(EXBTC); ...
diff(CHKTC); diff(DIATC); diff(GRATC); diff(NADTC); diff(missing_Tc)];
diff(CHKTC); diff(DIATC); diff(GRATC); diff(NADTC); diff(MISTC)];
TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/N_T,N_T]);
TIME_PER_STEP = sum(TIME_PER_FCT,2);
......@@ -57,7 +53,7 @@ checkfield_Ta = mean(diff(CHKTC));
grad_Ta = mean(diff(GRATC));
nadiab_Ta = mean(diff(NADTC));
diag_Ta = mean(diff(DIATC));
miss_Ta = mean(diff(missing_Tc));
miss_Ta = mean(diff(MISTC));
total_Ta = mean(diff(total_Tc));
names = {...
'Mrhs';
......@@ -79,9 +75,10 @@ Ts_A = [rhs_Ta adv_field_Ta ghost_Ta clos_Ta coll_Ta poisson_Ta...
NSTEP_PER_SAMP= mean(diff(TS0TC))/DTSIM;
%% Plots
fig = figure;
if 1
% subplot(121)
%% Area plot
fig = figure;
% colors = rand(N_T,3);
% colors = lines(N_T);
colors = distinguishable_colors(N_T);
......@@ -100,21 +97,21 @@ for i = 2:numel(x_)
yy_(2*i-1,:) = y_(i,:)/(dx/DTSIM);
yy_(2*i ,:) = y_(i,:)/(dx/DTSIM);
end
p1 = area(xx_,yy_,'LineStyle','none');
p1 = area(xx_/DTSIM,yy_,'LineStyle','none');
for i = 1:N_T; p1(i).FaceColor = colors(i,:);
% LEGEND{i} = sprintf('%s t=%1.1e[s] (%0.1f %s)',names{i},Ts_A(i),Ts_A(i)/total_Tc(end)*100,'\%');
LEGEND{i} = [names{i},' $\hat t=$',sprintf('%1.1e[s] (%0.1f %s)',Ts_A(i)/NSTEP_PER_SAMP,Ts_A(i)/total_Ta*100,'\%')];
end;
end
legend(LEGEND,'Location','bestoutside');
% legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Process', 'Missing')
xlabel('Sim. Time [$\rho_s/c_s$]'); ylabel('Step Comp. Time [s]')
xlim([TS0TC(2),TS0TC(end)]);
ylim([0, 1.1*CPUTI/(TS0TC(end)/DTSIM)])
xlabel('Simulation step'); ylabel('Comp. time per step [s]')
xlim([TS0TC(2),TS0TC(end)]./DTSIM);
ylim([0, 1.5*total_Ta/(dx/DTSIM)])
h_ = floor(CPUTI/3600);
m_ = floor(floor(CPUTI/60)-60*h_);
s_ = CPUTI - 3600*h_ - 60*m_;
title(sprintf('Gyacomo 2 (%.0f [h] ~%.0f [min] ~%.0f [s])',...
h_,m_,s_))
title([DATADIR,sprintf(' (%.0f [h] ~%.0f [min] ~%.0f [s])',...
h_,m_,s_)])
hold on
FIGNAME = 'profiler';
% save_figure
......@@ -145,20 +142,21 @@ FIGNAME = 'profiler';
end
if 0
subplot(122)
%% Histograms
fig = figure;
histogram(diff(rhs_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(adv_field_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(ghost_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(clos_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(coll_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(poisson_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(Sapj_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(process_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(checkfield_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(diag_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
% fig = figure;
histogram(diff(RHSTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(ADVTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(GHOTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(CLOTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(COLTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(POITC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(SAPTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(CHKTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(DIATC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
histogram(diff(MISTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
grid on;
legend('Compute RHS','Adv. fields','Ghosts comm', 'closure', 'collision','Poisson','Nonlin','Process','Check+sym', 'Diagnos.', 'Missing')
legend('Compute RHS','Adv. fields','Ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Missing')
xlabel('Step Comp. Time [s]'); ylabel('')
set(gca,'Xscale','log')
FIGNAME = 'profiler';
......
wk/HEL_lin_scan.m
+ 13
− 3
View file @ c36a0067
Nsim = 32;
gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory
addpath(genpath([gyacomodir,'matlab'])) % ... add
addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add
addpath(genpath([gyacomodir,'wk/parameters'])) % ... add
Nsim = 48;
gr_ = zeros(61,Nsim);
i = 1;
tau_a = logspace(0,-3,Nsim);
tau_a = logspace(0,-3.5,Nsim);
for i = 1:Nsim
run lin_Ivanov;
TAU = tau_a(i);
NU = 0.1*3/2/TAU/4; % Collision frequency
K_Ti = 0.35*2/TAU; % ion Temperature
DT = 1e-2;
DT = 5e-3;
PMAX = 2; JMAX = 1;
% PMAX = 2; JMAX = 0;
% PMAX = 0; JMAX = 1;
run fast_run;
% growth rates
[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
......
wk/fast_run.m 0 → 100644
+ 68
− 0
View file @ c36a0067
%% QUICK RUN SCRIPT
% This script runs a simulation directly from the Matlab framework.
% It is meant to run only small problems in linear for benchmarking and
% debugging purposes since it makes Matlab "busy".
%% Set up the paths for the necessary Matlab modules
gyacomodir = pwd;
gyacomodir = gyacomodir(1:end-2);
% mpirun = 'mpirun';
mpirun = '/opt/homebrew/bin/mpirun'; % for macos
addpath(genpath([gyacomodir,'matlab'])) % Add matlab folder
addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot folder
addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute folder
addpath(genpath([gyacomodir,'matlab/load'])) % Add load folder
addpath(genpath([gyacomodir,'wk/parameters'])) % Add parameters folder
%% Setup run or load an executable
RUN = 1; % To run or just to load
default_plots_options
% EXECNAME = 'gyacomo23_sp_save'; % single precision
% EXECNAME = 'gyacomo23_sp'; % single precision
EXECNAME = 'gyacomo23_dp'; % double precision
% EXECNAME = 'gyacomo23_debug'; % double precision
%% RUN
setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
RUN =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
% RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
% RUN =['time ',mpirun,' -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0;'];
% RUN =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;'];
% RUN =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
% RUN = ['./../../../bin/gyacomo23_sp 0;'];
MVOUT='cd ../../../wk;';
system([MVIN,RUN,MVOUT]);
end
%% quick load of basic data
filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
% Load outputs from jobnummin up to jobnummax
J0 = 0; J1 = 0;
data = {}; % Initialize data as an empty cell array
% load grids, inputs, and time traces
data = compile_results_low_mem(data,LOCALDIR,J0,J1);
if 0 % Activate or not
%% plot mode evolution and growth rates
% Load phi
[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
options.NORMALIZED = 0;
options.TIME = data.Ts3D;
% Time window to measure the growth of kx/ky modes
options.KY_TW = [0.5 1.0]*data.Ts3D(end);
options.KX_TW = [0.5 1.0]*data.Ts3D(end);
options.NMA = 1; % Set NMA option to 1
options.NMODES = 999; % Set how much modes we study
options.iz = 'avg'; % Compressing z
options.ik = 1; %
options.GOK2 = 0; % plot gamma/k^2
options.fftz.flag = 0; % Set fftz.flag option to 0
options.FIELD = 'phi';
options.SHOWFIG = 1;
[fig, wkykx, ekykx] = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
end
wk/get_miller_GENE_py.m
+ 1
− 1
View file @ c36a0067
......@@ -2,7 +2,7 @@ function [variables] = get_miller_GENE_py(prof_folder,rho)
filePath = [prof_folder,'/equilibrium.txt'];
command = ['python extract_miller_from_eqdsk.py ',...
command = ['python3 extract_miller_from_eqdsk.py ',...
filePath,' ',num2str(rho),' > tmp.txt'];
system(command);
......
wk/p3_geom_scan_analysis.m
+ 8
− 7
View file @ c36a0067
......@@ -14,7 +14,7 @@ GETFLUXSURFACE = 0;
% partition= '../results/paper_3/';
% Get the scan directory
switch 5
switch 6
case 1 % delta kappa scan
casename = 'DTT rho85';
partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/';
......@@ -56,20 +56,21 @@ switch 5
t1 = 50; t2 = 150;
case 5 % delta K_T tau=1
casename = 'DIIID rho95 $\tau=1$';
partition= '/misc/gyacomo23_outputs/paper_3/DIIID_tau_1_rho95_geom_scan/';
scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)';
partition= '../results/paper_3/DIIID_tau_1_rho95_geom_scan/';
% scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)';
scandir = '2_1_delta_RT_scan'; scanname= '(2,1)';
% scandir = '5x3x192x48x32_delta_RT_scan'; scanname= '(4,2)';
% scandir = 'delta_RT_scan_PJ_21'; scanname= '(2,1)';
nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1.0;
t1 = 200; t2 = 800;
t1 = 200; t2 = 500;
case 6 % delta K_T cold ions
casename = 'DIIID rho95 $\tau=10^{-3}$';
partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/';
partition= '../results/paper_3/DIIID_cold_ions_rho95_geom_scan/';
scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)';
nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1e3;
t1 = 200; t2 = 480;
nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1e3/2;
t1 = 150; t2 = 280;
case 7 % delta s_delta
casename = 'DIIID rho95 $\tau=10^{-3}$';
partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/';
......
wk/parameters/lin_Ivanov.m
+ 16
− 9
View file @ c36a0067
......@@ -3,12 +3,12 @@ SIMID = 'lin_Ivanov'; % Name of the simulation
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
TAU = 1e-3; % e/i temperature ratio
NU = 0.1/TAU; % Collision frequency
TAU = 0.001; % e/i temperature ratio
NU = 1.0*3/2/TAU/4; % Collision frequency
K_Ne = 0*2.22; % ele Density
K_Te = 0*6.96; % ele Temperature
K_Ni = 0*2.22; % ion Density gradient drive
K_Ti = 0.36*2/TAU; % ion Temperature
K_Ti = 1.0*2/TAU+0*5*TAU; % ion Temperature
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
NA = 1; % number of kinetic species
ADIAB_E = (NA==1); % adiabatic electron model
......@@ -20,9 +20,9 @@ J = 1;%P/2;
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
NX = 2; % real space x-gridpoints
NY = 40; % real space y-gridpoints
NY = 120; % real space y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
LY = 2*pi/0.15; % Size of the squared frequency domain in y direction
LY = 2*pi/0.05; % Size of the squared frequency domain in y direction
NZ = 1; % number of perpendicular planes (parallel grid)
SG = 0; % Staggered z grids option
NEXC = 0; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
......@@ -34,15 +34,18 @@ EPS = 0.0; % inverse aspect ratio
Q0 = 0.0; % safety factor
SHEAR = 0.0; % magnetic shear
KAPPA = 1.0; % elongation
S_KAPPA = 0.0;
DELTA = 0.0; % triangularity
S_DELTA = 0.0;
ZETA = 0.0; % squareness
S_ZETA = 0.0;
PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
SHIFT_Y = 0.0; % Shift in the periodic BC in z
NPOL = 1; % Number of poloidal turns
%% TIME PARAMETERS
TMAX = 50; % Maximal time unit
DT = 1e-2; % Time step
DT = 1e-3; % Time step
DTSAVE0D = 1; % Sampling time for 0D arrays
DTSAVE2D = -1; % Sampling time for 2D arrays
DTSAVE3D = 1; % Sampling time for 3D arrays
......@@ -51,12 +54,12 @@ JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
CO = 'SG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
GKCO = 0; % Gyrokinetic operator
CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
GKCO = 1; % Gyrokinetic operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; % Initialize zero-field quantities
% HRCY_CLOS = 'max_degree'; % Closure model for higher order moments
HRCY_CLOS = 'truncation'; % Closure model for higher order moments
HRCY_CLOS = 'max_degree'; % Closure model for higher order moments
DMAX = -1;
NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
NMAX = 0;
......@@ -93,3 +96,7 @@ BCKGD0 = 0.0; % Initial background
k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1.0; % Cutoff for collision operator
PB_PHASE = 0.0;
ADIAB_I = 0;
MHD_PD = 0;
EXBRATE = 0;
\ No newline at end of file
wk/plot_toroidal_fluxtube_geom.m
+ 6
− 6
View file @ c36a0067
......@@ -8,13 +8,13 @@ default_plots_options
Nx = 128;
Ny = 128;
Nz = 64;
Nturns = 1;
Nz = 128;
x = linspace(-60,60,Nx)*0.001;
y = linspace(-60,60,Ny)*0.001;
FIELD = ones(Nx,Ny,Nz);
z = linspace(-Nturns*pi,Nturns*pi,Nz);
N_field_lines = 120;
N_field_lines = 2;
N_magn_flux_surf = 1;
openangle = pi/3;
phi0 = openangle/2; phi1= 2*pi-openangle/2;
......@@ -28,9 +28,9 @@ select = 3;
switch select
case 1
% CBC
rho = 0.50; drho = 0.1;% Torus minor radius
rho = 0.3; drho = 0.1;% Torus minor radius
eps = 0.3;
q0 = 0.000001; % Flux tube safety factor
q0 = 1.4; % Flux tube safety factor
shear = 0.8;
kappa = 1.0;
s_kappa= 0.0;
......@@ -171,7 +171,7 @@ figure; set(gcf, 'Position', DIMENSIONS)
end
%plot field lines
if N_field_lines > 0
theta = linspace(-Nturns*pi, Nturns*pi, 256) ; % Poloidal angle
theta = linspace(-Nturns*pi, Nturns*pi, Nturns*256) ; % Poloidal angle
vecfl = Rvec(theta);
plot3(vecfl(1,:),vecfl(2,:),vecfl(3,:),'-b'); hold on;
% Multiple field lines
......@@ -183,7 +183,7 @@ figure; set(gcf, 'Position', DIMENSIONS)
y_ = [y_ sin(t)*vecfl(1,:) + cos(t)*vecfl(2,:)];
z_ = [z_ vecfl(3,:)];
end
plot3(x_,y_,z_,'.','color',[1.0 0.6 0.6]*0.8); hold on;
plot3(x_,y_,z_,'-','color',[1.0 0.6 0.6]*0.8); hold on;
end
%plot fluxe tube
if PLT_FTUBE
......
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