diff --git a/matlab/setup.m b/matlab/setup.m
index cf9e1f15734e23ee8cb8b8d20d997e8d5594220d..71b95cdc9136c34260cca173b6053437d6739e35 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -73,7 +73,7 @@ switch CO
end
COLL.coll_kcut = COLL_KCUT;
% Time integration and intialization parameters
-TIME_INTEGRATION.numerical_scheme = '''RK4''';
+TIME_INTEGRATION.numerical_scheme = ['''',NUMERICAL_SCHEME,''''];
INITIAL.INIT_OPT = ['''',INIT_OPT,''''];
INITIAL.ACT_ON_MODES = ['''',ACT_ON_MODES,''''];
INITIAL.init_background = BCKGD0;
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index e637d9dd64a137c983b4c1b7b083cd5039c7cd1b..b7bffafb4897af3121ce9c7601c18a8995eb61bd 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -64,7 +64,7 @@ options.INTERP = 1;
options.POLARPLOT = 0;
% options.NAME = '\phi';
% options.NAME = '\omega_z';
-options.NAME = 'N_e^{00}';
+options.NAME = 'N_i^{00}';
% options.NAME = 'v_x';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
@@ -95,7 +95,7 @@ options.NAME = '\phi';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
% options.NAME 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
@@ -189,8 +189,8 @@ end
if 0
%% Mode evolution
-options.NORMALIZED = 0;
-options.K2PLOT = [0.1 0.2 0.3 0.4];
+options.NORMALIZED = 1;
+options.K2PLOT = [0.1 0.5 1.0 2.0 6.0];
options.TIME = [00:1200];
options.NMA = 1;
options.NMODES = 5;
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 5e9a28c9d7d6d7dcfb740205eddcb8a67b992117..aa7ca49072a3e62716ca51e4e6a9eafde7ad28c7 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -62,10 +62,10 @@ gyacomodir = '/home/ahoffman/gyacomo/';
% resdir = 'testcases/miller_example';
% resdir = 'Miller/128x256x3x2_CBC_dt_5e-3';
%% CBC Miller
-% resdir = 'GCM_CBC/daint/Miller_GX_comparison';
+resdir = 'GCM_CBC/daint/Miller_GX_comparison';
%% RK3
-resdir = 'dbg/SSPx_RK3_test';
-resdir = 'dbg/SSPx_RK3_test/RK4';
-resdir = ['results/',resdir];
+% resdir = 'dbg/SSPx_RK3_test';
+% resdir = 'dbg/SSPx_RK3_test/RK4';
+% resdir = ['results/',resdir];
JOBNUMMIN = 00; JOBNUMMAX = 10;
run analysis_gyacomo
diff --git a/wk/lin_ETPY.m b/wk/lin_ETPY.m
index 9b6c383a95db4e8c877c6b1f79a8a3ceacc9e7c6..cdc2a5f583946c0b67b2ce17750125100d3135b7 100644
--- a/wk/lin_ETPY.m
+++ b/wk/lin_ETPY.m
@@ -3,57 +3,56 @@
% 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
+SIMID = 'lin_ETPY'; % 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_dbg';
-EXECNAME = 'gyacomo';
+EXECNAME = 'gyacomo_dbg';
+% EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.05; % Collision frequency
+NU = 0.01; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_Ne = 2.0; % ele Density '''
-K_Te = 0.4; % ele Temperature '''
-K_Ni = 2.0; % ion Density gradient drive
-K_Ti = 0.4; % ion Temperature '''
+K_Ne = 2.5; % ele Density '''
+K_Te = K_Ne/4; % ele Temperature '''
+K_Ni = K_Ne; % ion Density gradient drive
+K_Ti = K_Ni/4; % ion Temperature '''
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 = 8;
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 = 50; % '' y-gridpoints
+NY = 100; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
-LY = 2*pi/0.05; % Size of the squared frequency domain
+LY = 200;%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
%% GEOMETRY
-GEOMETRY= 'Z-pinch';
+GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
-KAPPA = 1.0; % elongation
+KAPPA = 0.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
%% TIME PARMETERS
-TMAX = 100; % Maximal time unit
+TMAX = 500; % 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
+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
SPSCP = 0; % Sampling per time unit for checkpoints
@@ -63,13 +62,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 = 0; % gyrokinetic operator
+GKCO = 1; % 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
+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
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
@@ -85,7 +85,7 @@ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
-MU_Z = 0.0; %
+MU_Z = 2.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index ebca810e6abee81ea9bba946a9a3cc6e9aa056b4..339b183351bc6791c41f482e36544a33835d264f 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -19,9 +19,9 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%% PHYSICAL PARAMETERS
NU = 0.05; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_Ne = 0*2.22; % ele Density '''
+K_Ne = 2.22; % ele Density '''
K_Te = 6.96; % ele Temperature '''
-K_Ni = 0*2.22; % ion Density gradient drive
+K_Ni = 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)
@@ -35,15 +35,15 @@ JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
NX = 20; % real space x-gridpoints
-NY = 2; % '' y-gridpoints
+NY = 20; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
-LY = 2*pi/0.3; % Size of the squared frequency domain
+LY = 2*pi/0.05; % Size of the squared frequency domain
NZ = 32; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
%% GEOMETRY
-% GEOMETRY= 's-alpha';
-GEOMETRY= 'miller';
+GEOMETRY= 's-alpha';
+% GEOMETRY= 'miller';
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
KAPPA = 1.0; % elongation
@@ -52,10 +52,10 @@ ZETA = 0.0; % squareness
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 25; % Maximal time unit
-DT = 3e-3; % Time step
+TMAX = 20; % 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
+SPS2D = -1; % Sampling per time unit for 2D arrays
SPS3D = 5; % Sampling per time unit for 2D arrays
SPS5D = 1/5; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
@@ -72,6 +72,7 @@ CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax)
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 = 'RK3'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
@@ -91,10 +92,11 @@ MU_Z = 2.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
-NOISE0 = 0.0e-5; % Init noise amplitude
-BCKGD0 = 1.0; % Init background
+NOISE0 = 1.0e-5; % Init noise amplitude
+BCKGD0 = 0.0; % Init background
GRADB = 1.0;
CURVB = 1.0;
+COLL_KCUT = 1000;
%%-------------------------------------------------------------------------
%% RUN
setup
diff --git a/wk/lin_RHT.m b/wk/lin_RHT.m
index dcdcc7529ab274daa0bd39c84c156e49604b2cf6..d9a230096da7fac82936e4e88166f68501fcd1fe 100644
--- a/wk/lin_RHT.m
+++ b/wk/lin_RHT.m
@@ -7,24 +7,24 @@ SIMID = 'lin_RHT'; % Name of the simulation
RUN = 1; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
-HELAZDIR = '/home/ahoffman/HeLaZ/';
-% EXECNAME = 'helaz3';
-EXECNAME = 'helaz3';
+GYACOMODIR = '/home/ahoffman/gyacomo/';
+% EXECNAME = 'gyacomo_dbg';
+EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.0; % Collision frequency
+NU = 0.1; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_Ne = 1.0; % ele Density '''
-K_Te = 4.0; % ele Temperature '''
-K_Ni = 2.0; % ion Density gradient drive
-K_Ti = 8.0; % ion Temperature '''
+K_Ne = 0.0; % ele Density '''
+K_Te = 0.0; % ele Temperature '''
+K_Ni = 0.0; % ion Density gradient drive
+K_Ti = 0.0; % ion Temperature '''
SIGMA_E = 0.5;%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
-BETA = 0.05; % electron plasma beta
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 0.00; % electron plasma beta
%% GRID PARAMETERS
P = 4;
J = P/2;
@@ -34,27 +34,29 @@ PMAXI = P; % " ions
JMAXI = J; % "
NX = 2; % real space x-gridpoints
NY = 2; % '' y-gridpoints
-LX = 2*pi/0.0628; % Size of the squared frequency domain
-LY = 2*pi/0.1; % Size of the squared frequency domain
+LX = pi/2.5; % Size of the squared frequency domain
+LY = 2*pi/0.5; % Size of the squared frequency domain
NZ = 16; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
%% GEOMETRY
% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
-GEOMETRY= 's-alpha';
-% GEOMETRY= 'circular';
+% GEOMETRY= 's-alpha';
+GEOMETRY= 'miller';
Q0 = 1.4; % safety factor
-SHEAR = 0.0; % magnetic shear
-NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+SHEAR = 0.8; % magnetic shear
+KAPPA = 1.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
+NEXC = 0; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) %0: adaptive
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 10; % Maximal time unit
-DT = 5e-3; % Time step
+TMAX = 50; % 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
-SPS3D = 10; % Sampling per time unit for 2D arrays
+SPS2D = -1; % Sampling per time unit for 2D arrays
+SPS3D = 5; % Sampling per time unit for 2D arrays
SPS5D = 1/5; % 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)
@@ -68,6 +70,7 @@ CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax)
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 = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
@@ -91,14 +94,14 @@ NOISE0 = 0.0e-5; % Init noise amplitude
BCKGD0 = 1.0; % Init background
GRADB = 1.0;
CURVB = 1.0;
+COLL_KCUT = 1000;
%%-------------------------------------------------------------------------
%% RUN
setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',GYACOMODIR,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
@@ -107,7 +110,7 @@ end
%% Load results
%%
filename = [SIMID,'/',PARAMS,'/'];
-LOCALDIR = [HELAZDIR,'results/',filename,'/'];
+LOCALDIR = [GYACOMODIR,'results/',filename,'/'];
% Load outputs from jobnummin up to jobnummax
JOBNUMMIN = 00; JOBNUMMAX = 00;
data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
@@ -164,13 +167,13 @@ end
if 0
%% Mode evolution
options.NORMALIZED = 0;
-options.K2PLOT = 1;
-options.TIME = [0:1000];
+options.K2PLOT = 2;
+options.TIME = [0:50];
options.NMA = 1;
options.NMODES = 1;
options.iz = 'avg';
fig = mode_growth_meter(data,options);
-save_figure(data,fig,'.png')
+% save_figure(data,fig,'.png')
end