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

parameters are in separate folders now

parent f6d580d1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 10 additions and 1020 deletions
wk/KBM_ky_PJ_scan.m deleted 100644 → 0
+ 0
− 218
View file @ f6d580d1
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';
EXECNAME = 'gyacomo23_sp';
% EXECNAME = 'gyacomo23_dp';
% EXECNAME = 'gyacomo23_debug';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
SIMID = 'lin_KBM'; % Name of the simulation
RERUN = 0; % rerun if the data does not exist
RUN = 1;
K_Ne = 3; % ele Density '''
K_Te = 4.5; % ele Temperature '''
K_Ni = 3; % ion Density gradient drive
K_Ti = 8; % ion Temperature '''
P_a = [2 4 8 16];
% P_a = 10;
ky_a= 0.05:0.1:0.75;
% ky_a = 0.05;
% collision setting
CO = 'DG';
GKCO = 1; % gyrokinetic operator
COLL_KCUT = 1.75;
NU = 1e-2;
% model
NA = 2;
BETA = 0.03; % electron plasma beta
% Geometry
GEOMETRY= 'miller';
% GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
DT0 = 5e-3;
TMAX = 15;
% arrays for the result
g_ky = zeros(numel(ky_a),numel(P_a),2);
g_avg= g_ky*0;
g_std= g_ky*0;
% Naming of the collision operator
if GKCO
CONAME = [CO,'GK'];
else
CONAME = [CO,'DK'];
end
j = 1;
for P = P_a
J = P/2;
i = 1;
for ky = ky_a
%% 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)
%% GRID PARAMETERS
DT = DT0/sqrt(J);
PMAX = P; % Hermite basis size
JMAX = P/2; % Laguerre "
NX = 12; % real space x-gridpoints
NY = 2;
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 2*pi/ky; % 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';
EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
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
DTSAVE0D = 1; % Sampling per time unit for 0D arrays
DTSAVE2D = -1; % Sampling per time unit for 2D arrays
DTSAVE3D = 0.5; % Sampling per time unit for 3D arrays
DTSAVE5D = 100; % Sampling per time unit for 5D arrays
JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% 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= '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 = 0;
W_GAMMA = 1; W_HF = 1;
W_PHI = 1; W_NA00 = 1;
W_DENS = 0; W_TEMP = 1;
W_NAPJ = 0; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
ADIAB_E = (NA==1);
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;
HYP_V = 'none';
HRCY_CLOS = 'truncation'; % Closure model for higher order moments
DMAX = -1;
NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
NMAX = 0;
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
k_gB = 1.0;
k_cB = 1.0;
%% RUN
setup
% naming
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
% check if data exist to run if no data
data_ = {};
try
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
Ntime = numel(data_.Ts0D);
catch
data_.outfilenames = [];
end
if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10)
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 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 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
end
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
if numel(data_.Ts3D)>10
if numel(data_.Ts3D)>5
% Load results after trying to run
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
% linear growth rate (adapted for 2D zpinch and fluxtube)
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
[~,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 = 2: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(it1:it2));
tw = double(data_.Ts3D(it1:it2));
% gr = polyfit(tw,to_measure,1);
gr = fit(tw,to_measure,'poly1');
err= confint(gr);
g_ky(i,j,ik) = gr.p1;
g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
end
[gmax, ikmax] = max(g_ky(i,j,:));
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
end
end
i = i + 1;
end
j = j + 1;
end
if 0
%% Check time evolution
figure;
to_measure = log(field_t);
plot(data_.Ts3D,to_measure); hold on
plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
end
%% take max growth rate among z coordinate
y_ = g_ky(:,:,2);
e_ = g_std(:,:,2);
%%
if(numel(ky_a)>1 && numel(P_a)>1)
%% Save metadata
pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
kymin = num2str(min(ky_a)); kymax= num2str(max(ky_a));
filename = [num2str(NX),'x',num2str(NZ),'_ky_',kymin,'_',kymax,...
'_P_',pmin,'_',pmax,'_',CONAME,'_',num2str(NU),'_be_',num2str(BETA),'.mat'];
metadata.name = filename;
metadata.kymin = ky;
metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),'$\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N)];
metadata.par = [num2str(NX),'x1x',num2str(NZ)];
metadata.nscan = 2;
metadata.s2name = '$P$';
metadata.s2 = P_a;
metadata.s1name = '$ky$';
metadata.s1 = ky_a;
metadata.dname = '$\gamma c_s/R$';
metadata.data = y_;
metadata.err = e_;
save([SIMDIR,filename],'-struct','metadata');
disp(['saved in ',SIMDIR,filename]);
clear metadata tosave
end
wk/KBM_load_metadata_scan.m deleted 100644 → 0
+ 0
− 172
View file @ f6d580d1
% Metadata path
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';
% datafname = 'lin_KBM/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.03.mat';
% datafname = 'lin_KBM/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.03.mat';
datafname ='lin_AUG_scan/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.000152.mat';
%% Chose if we filter gamma>0.05
FILTERGAMMA = 0;
%% Load data
fname = ['../results/',datafname];
d = load(fname);
if FILTERGAMMA
d.data = d.data.*(d.data>0.025);
d.err = d.err.*(d.data>0.025);
end
if 0
%% Pcolor of the peak
figure;
% [XX_,YY_] = meshgrid(d.s1,d.s2);
[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
pclr=imagesc_custom(XX_,YY_,d.data'.*(d.data>0)');
% pclr=contourf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
% pclr=surf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
title(d.title);
xlabel(d.s1name); ylabel(d.s2name);
set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
set(gca,'YTick',1:numel(d.s2),'YTicklabel',d.s2)
colormap(jet)
colormap(bluewhitered)
clb=colorbar;
clb.Label.String = '$\gamma c_s/R$';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
end
if 1
%% Scan along first dimension
figure
colors_ = jet(numel(d.s2));
for i = 1:numel(d.s2)
% plot(d.s1,d.data(:,i),'s-',...
plot(d.s1(d.data(:,i)>0),d.data((d.data(:,i)>0),i),'s-',...
'LineWidth',2.0,...
'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
'color',colors_(i,:));
% errorbar(d.s1,d.data(:,i),d.err(:,i),'s-',...
% 'LineWidth',2.0,...
% 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
% 'color',colors_(i,:));
hold on;
end
xlabel(d.s1name); ylabel(d.dname);title(d.title);
xlim([d.s1(1) d.s1(end)]);
colormap(colors_);
clb = colorbar;
% caxis([d.s2(1)-0.5,d.s2(end)+0.5]);
clim([1 numel(d.s2)+1]);
clb.Ticks=linspace(d.s2(1),d.s2(end),numel(d.s2));
clb.Ticks =1.5:numel(d.s2)+1.5;
clb.TickLabels=d.s2;
clb.Label.String = d.s2name;
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
end
if 0
%% Scan along second dimension
figure
colors_ = jet(numel(d.s1));
for i = 1:numel(d.s1)
plot(d.s2,d.data(i,:),'s-',...
'LineWidth',2.0,...
'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
'color',colors_(i,:));
% errorbar(d.s2,d.data(i,:),d.err(i,:),'s-',...
% 'LineWidth',2.0,...
% 'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
% 'color',colors_(i,:));
hold on;
end
xlabel(d.s2name); ylabel(d.dname);title(d.title);
xlim([d.s2(1) d.s2(end)]);
colormap(jet(numel(d.s1)));
clb = colorbar;
caxis([d.s1(1)-0.5,d.s1(end)+0.5]);
clb.Ticks=linspace(d.s1(1),d.s1(end),numel(d.s1));
clb.YTick=d.s1;
clb.Label.String = d.s1name;
clb.TickLabelInterpreter = 'latex';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
end
if 0
%% Convergence analysis
figure
% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
colors_ = jet(numel(d.s1));
for i = 1:numel(d.s1)
% target_ = d.data(i,end);
target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
if target_ > 0
eps_ = abs(target_ - d.data(i,1:end-1))/abs(target_);
semilogy(d.s2(1:end-1),eps_,'-s',...
'LineWidth',2.0,...
'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
'color',colors_(i,:));
hold on;
end
end
xlabel(d.s2name); ylabel('$\epsilon_r$');title(d.title);
xlim([d.s2(1) d.s2(end)]);
colormap(colors_);
clb = colorbar;
caxis([d.s1(1)-0.5,d.s1(end)+0.5]);
clb.Ticks=linspace(d.s1(1),d.s1(end),numel(d.s1));
clb.YTick=d.s1;
clb.Label.String = d.s1name;
clb.TickLabelInterpreter = 'latex';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
grid on;
end
if 0
%% Pcolor of the error
figure; i_ = 0;
% target_ = 2.72724991618068013377e-01; % Value for nuDGDK = 1.0, kT=6.96, (40,20), Nkx=8
% target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
% target_ = 2.71455e-01; % Value for nuDGDK = 0.001, kT=6.96, (40,20), Nkx=8
% target_ = 1.39427e-01; % Value for nuDGDK = 0.001, kT=5.3, (30,16), Nkx=8
target_ = 2.72510405826983714839e-01 % Value for nuDGDK = 0.001, kT=6.96, (50,25), Nkx=8
% target_ = 2.73048910051283844069e-01; % Value for nuDGDK = 0.0, kT=6.96, (40,20), Nkx=8
% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
% eps_ = log(abs(target_ - d.data)/abs(target_));
eps_ = max(-10,log(abs(target_ - d.data)/abs(target_)));
sign_ = 1;%sign(d.data - target_);
eps_ = d.data;
for i = 1:numel(d.s1)
for j = 1:numel(d.s2)
% target_ = d.data(i,end);
% target_ = d.data(i,end);
% target_ = d.data(end,j);
eps_(i,j) = log(abs(target_ - d.data(i,j))/target);
if target_ > 0
% eps_(i,:) = max(-12,log(abs(target_ - d.data(i,1:end))/abs(target_)));
% eps_(i,:) = log(abs(target_ - d.data(i,1:end))/abs(target_));
% eps_(i,:) = min(100,100*abs(target_ - d.data(i,1:end))/abs(target_));
else
end
end
end
[XX_,YY_] = meshgrid(d.s1,d.s2);
[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
pclr=imagesc_custom(XX_,YY_,eps_'.*(d.data>0)'.*sign_');
% pclr=contourf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_',5);
% pclr=surf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_');
title(d.title);
xlabel(d.s1name); ylabel(d.s2name);
set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
set(gca,'YTick',1:numel(d.s2),'YTicklabel',d.s2)
% colormap(jet)
colormap(bluewhitered)
% caxis([-10, 0]);
clb=colorbar;
clb.Label.String = '$\log(\epsilon_r)$';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
end
\ No newline at end of file
wk/lin_AUG.m deleted 100644 → 0
+ 0
− 183
View file @ f6d580d1
%% QUICK RUN SCRIPT
% This script creates a directory in /results and 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);
addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
%% Set simulation parameters
SIMID = 'lin_KBM'; % Name of the simulation
RUN = 1; % To run or just to load
default_plots_options
% EXECNAME = 'gyacomo23_sp'; % single precision
EXECNAME = 'gyacomo23_dp'; % double precision
% EXECNAME = 'gyacomo23_debug'; % single precision
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
NU = 0.1; % Collision frequency
TAU = 1.38; % e/i temperature ratio
K_Ne = 34; % ele Density '''
K_Te = 56; % ele Temperature '''
K_Ni = 34; % ion Density gradient drive
K_Ti = 21; % ion Temperature '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
NA = 2; % number of kinetic species
ADIAB_E = (NA==1); % adiabatic electron model
BETA = 1.52e-4; % electron plasma beta
MHD_PD = 0; % MHD pressure drift
NPOL = 1; % Number of poloidal turns
%% Set up grid parameters
P = 4;
J = P/2;%P/2;
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
NX = 6; % real space x-gridpoints
NY = 8; % real space y-gridpoints
LX = 2*pi/0.3; % Size of the squared frequency domain in x direction
LY = 2*pi/0.4; % Size of the squared frequency domain in y direction
NZ = 32*NPOL; % number of perpendicular planes (parallel grid)
SG = 0; % Staggered z grids option
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
%% GEOMETRY
% GEOMETRY= 's-alpha';
GEOMETRY= 'miller';
% GEOMETRY= 'circular';
EPS = 0.3; % inverse aspect ratio
Q0 = 5.7; % safety factor
SHEAR = 4.6;%*EPS; % magnetic shear
KAPPA = 1.8; % elongation
S_KAPPA = 0.0;
DELTA = 0.4; % 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
PB_PHASE= 0;
%% TIME PARAMETERS
TMAX = 5; % Maximal time unit
DT = 2e-3; % Time step
DTSAVE0D = 1; % Sampling per time unit for 0D arrays
DTSAVE2D = -1; % Sampling per time unit for 2D arrays
DTSAVE3D = 3; % Sampling per time unit for 3D arrays
DTSAVE5D = 100; % Sampling per time unit for 5D arrays
JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
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 = 'truncation'; % Closure model for higher order moments
DMAX = -1;
NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
NMAX = 0;
KERN = 0; % Kernel model (0 : GK)
INIT_OPT = 'mom00'; % Start simulation with a noisy mom00/phi/allmom
NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
%% OUTPUTS
W_DOUBLE = 1; % Output flag for double moments
W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
W_HF = 1; % Output flag for high-frequency potential energy
W_PHI = 1; % Output flag for potential
W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
W_DENS = 1; % Output flag for total density
W_TEMP = 1; % Output flag for temperature
W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% UNUSED PARAMETERS
% These parameters are usually not to play with in linear runs
MU = 0.1; % Hyperdiffusivity coefficient
MU_X = MU; % Hyperdiffusivity coefficient in x direction
MU_Y = MU; % Hyperdiffusivity coefficient in y direction
N_HD = 4; % Degree of spatial-hyperdiffusivity
MU_Z = 10.0; % Hyperdiffusivity coefficient in z direction
HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
LAMBDAD = 0.0; % Lambda Debye
NOISE0 = 0.0e-5; % Initial noise amplitude
BCKGD0 = 1.0e-5; % Initial background
k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1; % Cutoff for collision operator
ADIAB_I = 0; % adiabatic ion model
%%-------------------------------------------------------------------------
%% 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,' 1 6 1 0;'];
% RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
MVOUT='cd ../../../wk;';
system([MVIN,RUN,MVOUT]);
end
%% Analysis
% load
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.KX_TW = [0.2 1]*data.Ts3D(end);
options.KY_TW = [0.2 1]*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.fftz.flag = 0; % Set fftz.flag option to 0
fig = 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
if 1
%% Ballooning plot
[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
if data.inputs.BETA > 0
[data.PSI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'psi');
end
options.time_2_plot = [120];
options.kymodes = [0.1:0.1:1];
options.normalized = 1;
options.PLOT_KP = 0;
% options.field = 'phi';
fig = plot_ballooning(data,options);
title(GEOMETRY)
end
if 0
%% plot geom and metric coefficients
options.SHOW_FLUXSURF = 1;
options.SHOW_METRICS = 1;
[fig, geo_arrays] = plot_metric(data,options);
title(GEOMETRY)
end
\ No newline at end of file
wk/lin_DTT_AB_rho85_PT.m → wk/parameters/lin_DTT_AB_rho85_PT.m
+ 4
− 84
View file @ c2ec0ff1
%% QUICK RUN SCRIPT
% This script creates a directory in /results and 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';
addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
%% Set simulation parameters %% Set simulation parameters
SIMID = 'lin_DTT'; % Name of the simulation SIMID = 'lin_DTT_AB_rho85_PT'; % Name of the simulation
RUN = 1; % To run or just to load
default_plots_options
EXECNAME = 'gyacomo23_sp'; % single precision
% EXECNAME = 'gyacomo23_dp'; % double precision
%% Set up physical parameters %% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
NU = 0.05; % Collision frequency NU = 0.05; % Collision frequency
...@@ -38,11 +18,11 @@ P = 4; ...@@ -38,11 +18,11 @@ P = 4;
J = P/2;%P/2; J = P/2;%P/2;
PMAX = P; % Hermite basis size PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size JMAX = J; % Laguerre basis size
NX = 32; % real space x-gridpoints NX = 16; % real space x-gridpoints
NY = 16; % real space y-gridpoints NY = 16; % real space y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain in x direction LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
LY = 2*pi/0.2; % Size of the squared frequency domain in y direction LY = 2*pi/0.2; % Size of the squared frequency domain in y direction
NZ = 32; % number of perpendicular planes (parallel grid) NZ = 24; % number of perpendicular planes (parallel grid)
SG = 0; % Staggered z grids option SG = 0; % Staggered z grids option
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
...@@ -114,64 +94,4 @@ BCKGD0 = 0.0e-5; % Initial background ...@@ -114,64 +94,4 @@ BCKGD0 = 0.0e-5; % Initial background
k_gB = 1.0; % Magnetic gradient strength k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1; % Cutoff for collision operator COLL_KCUT = 1; % Cutoff for collision operator
ADIAB_I = 0; % adiabatic ion model ADIAB_I = 0; % adiabatic ion model
\ No newline at end of file
%%-------------------------------------------------------------------------
%% 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 8 ',gyacomodir,'bin/',EXECNAME,' 1 4 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
%% Analysis
% load
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.KX_TW = [0.2 1]*data.Ts3D(end);
options.KY_TW = [0.2 1]*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.fftz.flag = 0; % Set fftz.flag option to 0
fig = 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
if 1
%% Ballooning plot
[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
if data.inputs.BETA > 0
[data.PSI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'psi');
end
options.time_2_plot = [120];
options.kymodes = [0.25];
options.normalized = 1;
options.PLOT_KP = 0;
% options.field = 'phi';
fig = plot_ballooning(data,options);
end
wk/lin_ETG_adiab_i.m → wk/parameters/lin_ETG_adiab_i.m
+ 2
− 66
View file @ c2ec0ff1
%% QUICK RUN SCRIPT
% This script creates a directory in /results and 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);
addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
%% Set simulation parameters %% Set simulation parameters
SIMID = 'ETG_adiab_i'; % Name of the simulation SIMID = 'lin_ETG_adiab_i'; % Name of the simulation
RUN = 1; % To run or just to load
default_plots_options
EXECNAME = 'gyacomo23_sp'; % single precision
% EXECNAME = 'gyacomo23_dp'; % double precision
% EXECNAME = 'gyacomo23_debug'; % single precision
%% Set up physical parameters %% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
...@@ -113,50 +95,4 @@ NOISE0 = 1.0e-8; % Initial noise amplitude ...@@ -113,50 +95,4 @@ NOISE0 = 1.0e-8; % Initial noise amplitude
BCKGD0 = 0.0e-8; % Initial background BCKGD0 = 0.0e-8; % Initial background
k_gB = 1.0; % Magnetic gradient strength k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1; % Cutoff for collision operator COLL_KCUT = 1; % Cutoff for collision operator
\ No newline at end of file
%%-------------------------------------------------------------------------
%% 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,' 1 3 2 0;'];
% RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
MVOUT='cd ../../../wk;';
system([MVIN,RUN,MVOUT]);
end
%% Analysis
% load
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 1 % 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.KX_TW = [0.2 1]*data.Ts3D(end);
options.KY_TW = [0.2 1]*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.fftz.flag = 0; % Set fftz.flag option to 0
fig = 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/lin_Entropy.m → wk/parameters/lin_Entropy.m
+ 1
− 65
View file @ c2ec0ff1
%% QUICK RUN SCRIPT
% This script creates a directory in /results and 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);
addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
%% Set simulation parameters %% Set simulation parameters
SIMID = 'ETG_adiab_i'; % Name of the simulation SIMID = 'lin_Entropy'; % Name of the simulation
RUN = 1; % To run or just to load
default_plots_options
EXECNAME = 'gyacomo23_sp'; % single precision
% EXECNAME = 'gyacomo23_dp'; % double precision
% EXECNAME = 'gyacomo23_debug'; % single precision
%% Set up physical parameters %% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
...@@ -114,49 +96,3 @@ BCKGD0 = 0.0e-8; % Initial background ...@@ -114,49 +96,3 @@ BCKGD0 = 0.0e-8; % Initial background
k_gB = 1.0; % Magnetic gradient strength k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1; % Cutoff for collision operator COLL_KCUT = 1; % Cutoff for collision operator
%%-------------------------------------------------------------------------
%% 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 4 1 0;'];
% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 3 2 0;'];
% RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
MVOUT='cd ../../../wk;';
system([MVIN,RUN,MVOUT]);
end
%% Analysis
% load
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 1 % 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.KX_TW = [0.2 1]*data.Ts3D(end);
options.KY_TW = [0.2 1]*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.fftz.flag = 0; % Set fftz.flag option to 0
fig = 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/lin_ITG_salpha.m → wk/parameters/lin_ITG.m
+ 1
− 70
View file @ c2ec0ff1
%% QUICK RUN SCRIPT
% This script creates a directory in /results and 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);
addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
%% Set simulation parameters %% Set simulation parameters
SIMID = 'lin_ITG'; % Name of the simulation SIMID = 'lin_ITG'; % Name of the simulation
RUN = 1; % To run or just to load
default_plots_options
% EXECNAME = 'gyacomo23_dp'; % single precision
EXECNAME = 'gyacomo23_sp'; % double precision
% EXECNAME = 'gyacomo23_debug'; % double precision
%% Set up physical parameters %% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
...@@ -114,55 +96,4 @@ S_DELTA = 0.0; ...@@ -114,55 +96,4 @@ S_DELTA = 0.0;
S_ZETA = 0.0; S_ZETA = 0.0;
PB_PHASE= 0; PB_PHASE= 0;
ADIAB_I = 0; ADIAB_I = 0;
MHD_PD = 0; MHD_PD = 0;
%%------------------------------------------------------------------------- \ No newline at end of file
%% 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 =['mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
MVOUT='cd ../../../wk;';
system([MVIN,RUN,MVOUT]);
end
%% Analysis
% load
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
%% Plot heat flux evolution
figure
semilogy(data.Ts0D,data.HFLUX_X);
xlabel('$tc_s/R$'); ylabel('$Q_x$');
end
if 1 % 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.KX_TW = [0.2 1]*data.Ts3D(end);
options.KY_TW = [0.2 1]*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.fftz.flag = 0; % Set fftz.flag option to 0
fig = 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/lin_Ivanov.m → wk/parameters/lin_Ivanov.m
+ 2
− 87
View file @ c2ec0ff1
%% QUICK RUN SCRIPT
% This script creates a directory in /results and 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);
addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
%% Set simulation parameters %% Set simulation parameters
SIMID = 'dbg'; % Name of the simulation SIMID = 'lin_Ivanov'; % Name of the simulation
RUN = 1; % To run or just to load
default_plots_options
% EXECNAME = 'gyacomo23_dp'; % double precision
% To compile single precision gyacomo do 'make clean; make sp' in the /gyacomo folder
EXECNAME = 'gyacomo23_sp'; % single precision
% EXECNAME = 'gyacomo23_debug'; % single precision
%% Set up physical parameters %% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
TAU = 1e-3; % e/i temperature ratio TAU = 1e-3; % e/i temperature ratio
...@@ -111,70 +93,3 @@ BCKGD0 = 0.0; % Initial background ...@@ -111,70 +93,3 @@ BCKGD0 = 0.0; % Initial background
k_gB = 1.0; % Magnetic gradient strength k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1.0; % Cutoff for collision operator COLL_KCUT = 1.0; % Cutoff for collision operator
%%-------------------------------------------------------------------------
%% 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 4 1 0;'];
% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
% RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
MVOUT='cd ../../../wk;';
system([MVIN,RUN,MVOUT]);
end
%% Analysis
% load
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
%% Plot heat flux evolution
figure
semilogy(data.Ts0D,data.HFLUX_X);
xlabel('$tc_s/R$'); ylabel('$Q_x$');
end
if 1 % 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.KX_TW = [0.5 1]*data.Ts3D(end);
options.KY_TW = [0.5 1]*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.fftz.flag = 0; % Set fftz.flag option to 0
fig = 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
if 0
%% Hermite-Laguerre spectrum
[data.Napjz, data.Ts3D] = compile_results_3Da(LOCALDIR,J0,J1,'Napjz');
options.ST = 1;
options.NORMALIZED = 0;
options.LOGSCALE = 1;
options.FILTER = 0; %filter the 50% time-average of the spectrum from
fig = show_moments_spectrum(data,options);
end
% filename = '/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/Ivanov_2020_fig2_kT_0.26_chi_0.1.txt';
% % filename = '/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/Ivanov_2020_fig2_kT_1_chi_0.1.txt';
% dIV = load(filename);
%
% figure
% plot(dIV(:,1),2*dIV(:,2))
wk/lin_KBM.m → wk/parameters/lin_KBM.m
+ 0
− 75
View file @ c2ec0ff1
%% QUICK RUN SCRIPT
% This script creates a directory in /results and 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);
addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
%% Set simulation parameters %% Set simulation parameters
SIMID = 'lin_KBM'; % Name of the simulation SIMID = 'lin_KBM'; % Name of the simulation
RUN = 1; % To run or just to load
default_plots_options
EXECNAME = 'gyacomo23_sp'; % single precision
% EXECNAME = 'gyacomo23_dp'; % double precision
%% Set up physical parameters %% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
...@@ -108,61 +91,3 @@ BCKGD0 = 1.0e-5; % Initial background ...@@ -108,61 +91,3 @@ BCKGD0 = 1.0e-5; % Initial background
k_gB = 1.0; % Magnetic gradient strength k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1; % Cutoff for collision operator COLL_KCUT = 1; % Cutoff for collision operator
%%-------------------------------------------------------------------------
%% 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,' 1 6 1 0;'];
% RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
MVOUT='cd ../../../wk;';
system([MVIN,RUN,MVOUT]);
end
%% Analysis
% load
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.KX_TW = [0.2 1]*data.Ts3D(end);
options.KY_TW = [0.2 1]*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.fftz.flag = 0; % Set fftz.flag option to 0
fig = 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
if 1
%% Ballooning plot
[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
if data.inputs.BETA > 0
[data.PSI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'psi');
end
options.time_2_plot = [120];
options.kymodes = [0.25];
options.normalized = 1;
% options.field = 'phi';
fig = plot_ballooning(data,options);
end
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