diff --git a/matlab/MOLI_kperp_scan.m b/matlab/MOLI_kperp_scan.m
deleted file mode 100644
index 9920550ee6c1be9b32e97ad9ef36d1ff8f6e7167..0000000000000000000000000000000000000000
--- a/matlab/MOLI_kperp_scan.m
+++ /dev/null
@@ -1,166 +0,0 @@
-%% Run MOLI for a scan in kperp
-%% Move to MOLI workspace
-cd ../../MoliSolver/MOLI/workspace/
-%% Add paths
-my_paths
-
-%% Directories
-ROOT = '/home/ahoffman/Documents/MoliSolver';
-options.dirs.COSOlverdir = fullfile(ROOT,'COSOlver');
-options.dirs.MOLIdir = fullfile(ROOT,'MOLI');
-
-%% MOLI Physical and Main Parameters
-
-% Solve DK AND/OR GK Linear Moment Hierarchy.
-options.DKI = 0; % First-order DK
-options.DKII = 0; % Second-order DK
-options.GK = 1; % Gyrokinetic GK
-options.EM = 0; % Include Electromagnetic effects (only for GK=1)
-options.GD = 0; % Gyro-Drift
-options.GDI = 0; % First/second order
-
-% Solve MOLI
-options.MOLI = 1; % 1 -> Solve MOLI, 0 -> off
-
-% MOLI Solver
-options.solver.solver = 3;
-
-% MOLI Linear Fit Solver
-options.LinFitSolver = 2;
-
-%% Main parameter scan
-
-% Closure by truncation
-params.Pmaxi = GRID.pmaxi; % parallel ion Hermite moments
-params.Jmaxi = GRID.jmaxi; % perpendicular ion Laguerre moments
-params.Pmaxe = GRID.pmaxe; % parallel electron Hermite moments
-params.Jmaxe = GRID.jmaxe; % perpendicular electron Laguerre moments
-
-% w/wo ions
-options.ions = 1; % if ions are present -> 1, 0 otherwise
-
-% Adiabatic electrons
-options.electrons = 1; % 0 -> adiabatic electrons, 1 no adiabatic electrons
-
-% w/wo soundwaves
-options.sw = 1; % 1 -> sound waves on, 0 -> put ion parallel velocity row/column to 0
-
-%% Collision Operator Models and COSOlver Input Parameters
-options.collI = MODEL.CO; % collI=-2 -> Dougherty, -1 -> COSOlver, 0 -> Lenard-Bernstein, other -> hyperviscosity exponent
-options.collGK = 0; % collDKGK =1 -> GK collision operator, else DK collision operator
-options.COSOlver.GKE = 0;
-options.COSOlver.GKI = 0;
-
-% COSOlver Input Parameters (if collI = -1 only)
-options.COSOlver.eecolls = 1; % 1 -> electron-electron collisions, 0 -> off
-options.COSOlver.iicolls = 1; % 1 -> ion-ion collisions, 0 -> off
-options.COSOlver.eicolls = 1; % 1 -> electron-ion collisions (e-i) on, 0 -> off
-options.COSOlver.iecolls = 1; % 1 -> ion-electron collisions (i-e) on, 0 -> off
-
-% Collisional Coulomb sum bounds (only if collI = -1, i.e. Coulomb)
-options.COSOlver.lmaxx = 10; % upper bound collision operator first sum first species
-options.COSOlver.kmaxx = 10; % upper bound collision operator second sum first species
-options.COSOlver.nmaxx = options.COSOlver.lmaxx; % upper bound collision operator first sum second species
-options.COSOlver.qmaxx = options.COSOlver.kmaxx; % upper bound collision operator second sum second species
-
-% Collsion FLR sum bounds
-options.COSOlver.nemaxxFLR = 0; % upper bound FLR electron collision
-options.COSOlver.nimaxxFLR = 0; % upper bound FLR ion collision
-
-% Collision Operator Model
-% Set electron/ion test and back-reaction model operator
-%
-% 0 => Coulomb Collisions
-options.COSOlver.ETEST = 1; % 0 --> Buffer Operator, 1 --> Coulomb, 2 --> Lorentz
-options.COSOlver.EBACK = 1;
-options.COSOlver.ITEST = 1;
-options.COSOlver.IBACK = 1;
-options.COSOlver.ESELF = 1;
-options.COSOlver.ISELF = 1;
-
-options.COSOlver.OVERWRITE = 0; % overwrite collisional matrices even if exists
-
-options.COSOlver.cmd = 'mpirun -np 6 ./bin/CO 2 2 2';
-%% Physical Parameters
-
-% Toroidal effects
-options.magnetic = 1; % 1-> Add toroidal magnetic gradient drift resonance effects
-
-% Physical Parameters
-params.tau = MODEL.tau_i; % Ti/Te
-params.nu = MODEL.nu; % electron/ion collision frequency ... only for nu/ omega_pe < nuoveromegapemax (electron plasma frequency) [See Banks et al. (2017)]
-params.nuoveromegapemax = inf; % Maximum ratio between electron/ion collision frequency and electron plasma frequency [See Banks et al. (2017)]. Set to inf if not desired !!!
-params.mu = MODEL.sigma_e; % sqrt(m_e/m_i)
-params.kpar = 0.0; % normalized parallel wave number to the major radius
-params.kperp = GRID.kzmin; % normalized perpendicular toroidal wave number to the soundLarmor radius. Note: If ions ==0 (e.g. EPW), kperp --> b
-params.kr = GRID.krmin; % Radial component of perpendicular vector
-params.alphaD = 0.0; % (k*Debye length)^2
-params.Rn = MODEL.eta_n; % Major Radius / Background density gradient length
-params.RTe = MODEL.eta_T; % Major Radius * normalized kperp / Background electron temperature gradient length
-params.RTi = MODEL.eta_T; % Major Radius * normalized kperp / Background ion temperature gradient length
-params.Rphi = 0.0; % Major Radius * normalized kperp / Background potentiel gradient length [presence of shear] - only for GK
-params.betae = 1e-6; % Electron Beta plasma.
-
-params.rhostar = 1e-5; % sound Larmor Radius/Major Radius ~ sqrt(Te)/(R_0*B).
-params.n0 = INITIAL.init_background; % initial density perturbation
-
-params.gradB = MODEL.eta_B; % Magnetic field gradient
-params.curvB = MODEL.eta_B; % Curvature of B
-params.trappB = 0.0; % Trapping term
-
-%% MOLI Options
-
-% Save data in results dir
-options.save = 0;
-options.verbose = 0;
-options.dbg = 0;
-
-options.DR = 0; % 1 -> Solve kinetic dispersion relation,
-options.KineticDR = 0; % Solve kinetic dispersion relation (Landau integral) for the given theory
-
-% Compute the kinetic susceptibility for EPW only
-options.SPTBLTY = 0;
-
-options.nharm = 1; % Number of harmonics in disp. rel. 1 and 4
-
-wlim = 5.0;
-options.DRsolver.wr_min = -wlim; % Minimum of real part.
-options.DRsolver.wr_max = wlim; % Maximum of real part.
-options.DRsolver.wi_min = -wlim; % Minimum of imag part.
-options.DRsolver.wi_max = wlim; % Maximum of imag part.
-options.DRsolver.nw = 300; % Grid resolution
-
-% Disp. Rel. Options
-options.FLRmodel = 0; % 1 -> Truncated Laguerre, 0 -> Exact representation
-options.FluidLandau = 0; % 1 -> Add Landau Fluid Closure to Fluid Dispersion Relation, 0 -> off
-options.deltaLandau = 0; % 1 -> Hammet-Perkins closure on, 0 -> off
-
-% Fluid dispersion relation
-options.FluidDR = 0; % Solve annamaria's fluid equations
-options.Fluid.sITGmm = 0;
-
-% Define scan parameters
-options.fscan = 1; % 1 -> peform scan over scan.list, 0-> off
-
-options.scan.list = {'kperp'};% List of scan parameters. If empty, solve MOLI with params
-
-options.scan.kperp.array = linspace(GRID.kzmin, GRID.kzmax, GRID.nkz);
-
-% Time-Evolution Problem [Solver==3] ...
-options.solver.TimeSolver.dt = BASIC.dt; % timestep of time evolution (R/c_s or 1/(k v/the) units)
-options.solver.TimeSolver.tmax = BASIC.tmax;
-options.solver.TimeSolver.Trun = BASIC.tmax; % total time to run time evolution
-options.solver.TimeSolver.t_fit_min = 0.05; % Phase-Mixing fit Lower time limit
-options.solver.TimeSolver.t_fit_max = 8; % Phase-Mixing fit Upper time limit
-options.solver.TimeSolver.en_fit_min = 0.15; % Entropy Mode fit Lower time limit
-options.solver.TimeSolver.en_fit_max = 0.3; % Entropy Mode fit Upper time limit
-options.solver.TimeSolver.movie = 0; % Display movie if 1, last frame otherwise
-options.solver.TimeSolver.save = 0; % 1 --> save during fscan, Warning: memory storage
-
-
-%% Run MOLI
-% Solve the MOLI
-[results,params,options] = MOLI_Control(params,options);
-
-%% Return to HeLaZ workspace
-cd ../../../HeLaZ/wk
diff --git a/matlab/MOLI_time_solver.m b/matlab/MOLI_time_solver.m
deleted file mode 100644
index c6c6ed309ca332bff8c79250e30b3f22df865268..0000000000000000000000000000000000000000
--- a/matlab/MOLI_time_solver.m
+++ /dev/null
@@ -1,164 +0,0 @@
-%% Run MOLI for a time evolution of the moments at a given kperp
-%% Move to MOLI workspace
-cd ../../MoliSolver/MOLI/workspace/
-%% Add paths
-my_paths
-
-%% Directories
-ROOT = '/home/ahoffman/Documents/MoliSolver';
-options.dirs.COSOlverdir = fullfile(ROOT,'COSOlver');
-options.dirs.MOLIdir = fullfile(ROOT,'MOLI');
-
-%% MOLI Physical and Main Parameters
-
-% Solve DK AND/OR GK Linear Moment Hierarchy.
-options.DKI = 0; % First-order DK
-options.DKII = 0; % Second-order DK
-options.GK = 1; % Gyrokinetic GK
-options.EM = 0; % Include Electromagnetic effects (only for GK=1)
-options.GD = 0; % Gyro-Drift
-options.GDI = 0; % First/second order
-
-% Solve MOLI
-options.MOLI = 1; % 1 -> Solve MOLI, 0 -> off
-
-% MOLI Solver
-options.solver.solver = 3;
-
-% MOLI Linear Fit Solver
-options.LinFitSolver = 0;
-
-%% Main parameter scan
-
-% Closure by truncation
-params.Pmaxi = GRID.pmaxi; % parallel ion Hermite moments
-params.Jmaxi = GRID.jmaxi; % perpendicular ion Laguerre moments
-params.Pmaxe = GRID.pmaxe; % parallel electron Hermite moments
-params.Jmaxe = GRID.jmaxe; % perpendicular electron Laguerre moments
-
-% w/wo ions
-options.ions = 1; % if ions are present -> 1, 0 otherwise
-
-% Adiabatic electrons
-options.electrons = 1; % 0 -> adiabatic electrons, 1 no adiabatic electrons
-
-% w/wo soundwaves
-options.sw = 1; % 1 -> sound waves on, 0 -> put ion parallel velocity row/column to 0
-
-%% Collision Operator Models and COSOlver Input Parameters
-options.collI = MODEL.CO; % collI=-2 -> Dougherty, -1 -> COSOlver, 0 -> Lenard-Bernstein, other -> hyperviscosity exponent
-options.collGK = 0; % collDKGK =1 -> GK collision operator, else DK collision operator
-options.COSOlver.GKE = 0;
-options.COSOlver.GKI = 0;
-
-% COSOlver Input Parameters (if collI = -1 only)
-options.COSOlver.eecolls = 1; % 1 -> electron-electron collisions, 0 -> off
-options.COSOlver.iicolls = 1; % 1 -> ion-ion collisions, 0 -> off
-options.COSOlver.eicolls = 1; % 1 -> electron-ion collisions (e-i) on, 0 -> off
-options.COSOlver.iecolls = 1; % 1 -> ion-electron collisions (i-e) on, 0 -> off
-
-% Collisional Coulomb sum bounds (only if collI = -1, i.e. Coulomb)
-options.COSOlver.lmaxx = 10; % upper bound collision operator first sum first species
-options.COSOlver.kmaxx = 10; % upper bound collision operator second sum first species
-options.COSOlver.nmaxx = options.COSOlver.lmaxx; % upper bound collision operator first sum second species
-options.COSOlver.qmaxx = options.COSOlver.kmaxx; % upper bound collision operator second sum second species
-
-% Collsion FLR sum bounds
-options.COSOlver.nemaxxFLR = 0; % upper bound FLR electron collision
-options.COSOlver.nimaxxFLR = 0; % upper bound FLR ion collision
-
-% Collision Operator Model
-% Set electron/ion test and back-reaction model operator
-%
-% 0 => Coulomb Collisions
-options.COSOlver.ETEST = 1; % 0 --> Buffer Operator, 1 --> Coulomb, 2 --> Lorentz
-options.COSOlver.EBACK = 1;
-options.COSOlver.ITEST = 1;
-options.COSOlver.IBACK = 1;
-options.COSOlver.ESELF = 1;
-options.COSOlver.ISELF = 1;
-
-options.COSOlver.OVERWRITE = 0; % overwrite collisional matrices even if exists
-
-options.COSOlver.cmd = 'mpirun -np 6 ./bin/CO 2 2 2';
-%% Physical Parameters
-
-% Toroidal effects
-options.magnetic = 1; % 1-> Add toroidal magnetic gradient drift resonance effects
-
-% Physical Parameters
-params.tau = MODEL.tau_i; % Ti/Te
-params.nu = MODEL.nu; % electron/ion collision frequency ... only for nu/ omega_pe < nuoveromegapemax (electron plasma frequency) [See Banks et al. (2017)]
-params.nuoveromegapemax = inf; % Maximum ratio between electron/ion collision frequency and electron plasma frequency [See Banks et al. (2017)]. Set to inf if not desired !!!
-params.mu = MODEL.sigma_e; % sqrt(m_e/m_i)
-params.kpar = 0.0; % normalized parallel wave number to the major radius
-params.kperp = kz_MOLI; % normalized perpendicular toroidal wave number to the soundLarmor radius. Note: If ions ==0 (e.g. EPW), kperp --> b
-params.kr = kr_MOLI; % Radial component of perpendicular vector
-params.alphaD = 0.0; % (k*Debye length)^2
-params.Rn = MODEL.eta_n; % Major Radius / Background density gradient length
-params.RTe = MODEL.eta_T; % Major Radius * normalized kperp / Background electron temperature gradient length
-params.RTi = MODEL.eta_T; % Major Radius * normalized kperp / Background ion temperature gradient length
-params.Rphi = 0.0; % Major Radius * normalized kperp / Background potentiel gradient length [presence of shear] - only for GK
-params.betae = 1e-6; % Electron Beta plasma.
-
-params.rhostar = 1e-5; % sound Larmor Radius/Major Radius ~ sqrt(Te)/(R_0*B).
-params.n0 = INITIAL.init_background; % initial density perturbation
-
-params.gradB = MODEL.eta_B; % Magnetic field gradient
-params.curvB = MODEL.eta_B; % Curvature of B
-params.trappB = 0.0; % Trapping term
-
-%% MOLI Options
-
-% Save data in results dir
-options.save = 0;
-options.verbose = 0;
-options.dbg = 0;
-
-options.DR = 0; % 1 -> Solve kinetic dispersion relation,
-options.KineticDR = 0; % Solve kinetic dispersion relation (Landau integral) for the given theory
-
-% Compute the kinetic susceptibility for EPW only
-options.SPTBLTY = 0;
-
-options.nharm = 1; % Number of harmonics in disp. rel. 1 and 4
-
-wlim = 5.0;
-options.DRsolver.wr_min = -wlim; % Minimum of real part.
-options.DRsolver.wr_max = wlim; % Maximum of real part.
-options.DRsolver.wi_min = -wlim; % Minimum of imag part.
-options.DRsolver.wi_max = wlim; % Maximum of imag part.
-options.DRsolver.nw = 300; % Grid resolution
-
-% Disp. Rel. Options
-options.FLRmodel = 0; % 1 -> Truncated Laguerre, 0 -> Exact representation
-options.FluidLandau = 0; % 1 -> Add Landau Fluid Closure to Fluid Dispersion Relation, 0 -> off
-options.deltaLandau = 0; % 1 -> Hammet-Perkins closure on, 0 -> off
-
-% Fluid dispersion relation
-options.FluidDR = 0; % Solve annamaria's fluid equations
-options.Fluid.sITGmm = 0;
-
-% Define scan parameters
-options.fscan = 0; % 1 -> peform scan over scan.list, 0-> off
-
-options.scan.list = {};% List of scan parameters. If empty, solve MOLI with params
-
-% Time-Evolution Problem [Solver==3] ...
-options.solver.TimeSolver.dt = BASIC.dt; % timestep of time evolution (R/c_s or 1/(k v/the) units)
-options.solver.TimeSolver.tmax = BASIC.tmax;
-options.solver.TimeSolver.Trun = BASIC.tmax; % total time to run time evolution
-options.solver.TimeSolver.t_fit_min = 0.05; % Phase-Mixing fit Lower time limit
-options.solver.TimeSolver.t_fit_max = 8; % Phase-Mixing fit Upper time limit
-options.solver.TimeSolver.en_fit_min = 0.15; % Entropy Mode fit Lower time limit
-options.solver.TimeSolver.en_fit_max = 0.3; % Entropy Mode fit Upper time limit
-options.solver.TimeSolver.movie = 0; % Display movie if 1, last frame otherwise
-options.solver.TimeSolver.save = 0; % 1 --> save during fscan, Warning: memory storage
-
-
-%% Run MOLI
-% Solve the MOLI
-[results,params,options] = MOLI_Control(params,options);
-
-%% Return to HeLaZ workspace
-cd ../../../HeLaZ/wk
diff --git a/matlab/MOLI_time_solver_2D.m b/matlab/MOLI_time_solver_2D.m
deleted file mode 100644
index a356c5e71f300c24d201f755985c6c50217c2ec1..0000000000000000000000000000000000000000
--- a/matlab/MOLI_time_solver_2D.m
+++ /dev/null
@@ -1,164 +0,0 @@
-%% Run MOLI for a time evolution of the moments at a given kperp
-%% Move to MOLI workspace
-cd ../../MoliSolver/MOLI/workspace/
-%% Add paths
-my_paths
-
-%% Directories
-ROOT = '/home/ahoffman/Documents/MoliSolver';
-options.dirs.COSOlverdir = fullfile(ROOT,'COSOlver');
-options.dirs.MOLIdir = fullfile(ROOT,'MOLI');
-
-%% MOLI Physical and Main Parameters
-
-% Solve DK AND/OR GK Linear Moment Hierarchy.
-options.DKI = 0; % First-order DK
-options.DKII = 0; % Second-order DK
-options.GK = 1; % Gyrokinetic GK
-options.EM = 0; % Include Electromagnetic effects (only for GK=1)
-options.GD = 0; % Gyro-Drift
-options.GDI = 0; % First/second order
-
-% Solve MOLI
-options.MOLI = 1; % 1 -> Solve MOLI, 0 -> off
-
-% MOLI Solver
-options.solver.solver = 3;
-
-% MOLI Linear Fit Solver
-options.LinFitSolver = 0;
-
-%% Main parameter scan
-
-% Closure by truncation
-params.Pmaxi = GRID.pmaxi; % parallel ion Hermite moments
-params.Jmaxi = GRID.jmaxi; % perpendicular ion Laguerre moments
-params.Pmaxe = GRID.pmaxe; % parallel electron Hermite moments
-params.Jmaxe = GRID.jmaxe; % perpendicular electron Laguerre moments
-
-% w/wo ions
-options.ions = 1; % if ions are present -> 1, 0 otherwise
-
-% Adiabatic electrons
-options.electrons = 1; % 0 -> adiabatic electrons, 1 no adiabatic electrons
-
-% w/wo soundwaves
-options.sw = 1; % 1 -> sound waves on, 0 -> put ion parallel velocity row/column to 0
-
-%% Collision Operator Models and COSOlver Input Parameters
-options.collI = MODEL.CO; % collI=-2 -> Dougherty, -1 -> COSOlver, 0 -> Lenard-Bernstein, other -> hyperviscosity exponent
-options.collGK = 0; % collDKGK =1 -> GK collision operator, else DK collision operator
-options.COSOlver.GKE = 0;
-options.COSOlver.GKI = 0;
-
-% COSOlver Input Parameters (if collI = -1 only)
-options.COSOlver.eecolls = 1; % 1 -> electron-electron collisions, 0 -> off
-options.COSOlver.iicolls = 1; % 1 -> ion-ion collisions, 0 -> off
-options.COSOlver.eicolls = 1; % 1 -> electron-ion collisions (e-i) on, 0 -> off
-options.COSOlver.iecolls = 1; % 1 -> ion-electron collisions (i-e) on, 0 -> off
-
-% Collisional Coulomb sum bounds (only if collI = -1, i.e. Coulomb)
-options.COSOlver.lmaxx = 10; % upper bound collision operator first sum first species
-options.COSOlver.kmaxx = 10; % upper bound collision operator second sum first species
-options.COSOlver.nmaxx = options.COSOlver.lmaxx; % upper bound collision operator first sum second species
-options.COSOlver.qmaxx = options.COSOlver.kmaxx; % upper bound collision operator second sum second species
-
-% Collsion FLR sum bounds
-options.COSOlver.nemaxxFLR = 0; % upper bound FLR electron collision
-options.COSOlver.nimaxxFLR = 0; % upper bound FLR ion collision
-
-% Collision Operator Model
-% Set electron/ion test and back-reaction model operator
-%
-% 0 => Coulomb Collisions
-options.COSOlver.ETEST = 1; % 0 --> Buffer Operator, 1 --> Coulomb, 2 --> Lorentz
-options.COSOlver.EBACK = 1;
-options.COSOlver.ITEST = 1;
-options.COSOlver.IBACK = 1;
-options.COSOlver.ESELF = 1;
-options.COSOlver.ISELF = 1;
-
-options.COSOlver.OVERWRITE = 0; % overwrite collisional matrices even if exists
-
-options.COSOlver.cmd = 'mpirun -np 6 ./bin/CO 2 2 2';
-%% Physical Parameters
-
-% Toroidal effects
-options.magnetic = 1; % 1-> Add toroidal magnetic gradient drift resonance effects
-
-% Physical Parameters
-params.tau = MODEL.tau_i; % Ti/Te
-params.nu = MODEL.nu; % electron/ion collision frequency ... only for nu/ omega_pe < nuoveromegapemax (electron plasma frequency) [See Banks et al. (2017)]
-params.nuoveromegapemax = inf; % Maximum ratio between electron/ion collision frequency and electron plasma frequency [See Banks et al. (2017)]. Set to inf if not desired !!!
-params.mu = MODEL.sigma_e; % sqrt(m_e/m_i)
-params.kpar = 0.0; % normalized parallel wave number to the major radius
-params.kperp = kz; % normalized perpendicular toroidal wave number to the soundLarmor radius. Note: If ions ==0 (e.g. EPW), kperp --> b
-params.kr = kr; % Radial component of perpendicular vector
-params.alphaD = 0.0; % (k*Debye length)^2
-params.Rn = MODEL.eta_n; % Major Radius / Background density gradient length
-params.RTe = MODEL.eta_T; % Major Radius * normalized kperp / Background electron temperature gradient length
-params.RTi = MODEL.eta_T; % Major Radius * normalized kperp / Background ion temperature gradient length
-params.Rphi = 0.0; % Major Radius * normalized kperp / Background potentiel gradient length [presence of shear] - only for GK
-params.betae = 1e-6; % Electron Beta plasma.
-
-params.rhostar = 1e-5; % sound Larmor Radius/Major Radius ~ sqrt(Te)/(R_0*B).
-params.n0 = INITIAL.init_background; % initial density perturbation
-
-params.gradB = MODEL.eta_B; % Magnetic field gradient
-params.curvB = MODEL.eta_B; % Curvature of B
-params.trappB = 0.0; % Trapping term
-
-%% MOLI Options
-
-% Save data in results dir
-options.save = 0;
-options.verbose = 0;
-options.dbg = 0;
-
-options.DR = 0; % 1 -> Solve kinetic dispersion relation,
-options.KineticDR = 0; % Solve kinetic dispersion relation (Landau integral) for the given theory
-
-% Compute the kinetic susceptibility for EPW only
-options.SPTBLTY = 0;
-
-options.nharm = 1; % Number of harmonics in disp. rel. 1 and 4
-
-wlim = 5.0;
-options.DRsolver.wr_min = -wlim; % Minimum of real part.
-options.DRsolver.wr_max = wlim; % Maximum of real part.
-options.DRsolver.wi_min = -wlim; % Minimum of imag part.
-options.DRsolver.wi_max = wlim; % Maximum of imag part.
-options.DRsolver.nw = 300; % Grid resolution
-
-% Disp. Rel. Options
-options.FLRmodel = 0; % 1 -> Truncated Laguerre, 0 -> Exact representation
-options.FluidLandau = 0; % 1 -> Add Landau Fluid Closure to Fluid Dispersion Relation, 0 -> off
-options.deltaLandau = 0; % 1 -> Hammet-Perkins closure on, 0 -> off
-
-% Fluid dispersion relation
-options.FluidDR = 0; % Solve annamaria's fluid equations
-options.Fluid.sITGmm = 0;
-
-% Define scan parameters
-options.fscan = 0; % 1 -> peform scan over scan.list, 0-> off
-
-options.scan.list = {};% List of scan parameters. If empty, solve MOLI with params
-
-% Time-Evolution Problem [Solver==3] ...
-options.solver.TimeSolver.dt = BASIC.dt; % timestep of time evolution (R/c_s or 1/(k v/the) units)
-options.solver.TimeSolver.tmax = BASIC.tmax;
-options.solver.TimeSolver.Trun = BASIC.tmax; % total time to run time evolution
-options.solver.TimeSolver.t_fit_min = 0.05; % Phase-Mixing fit Lower time limit
-options.solver.TimeSolver.t_fit_max = 8; % Phase-Mixing fit Upper time limit
-options.solver.TimeSolver.en_fit_min = 0.15; % Entropy Mode fit Lower time limit
-options.solver.TimeSolver.en_fit_max = 0.3; % Entropy Mode fit Upper time limit
-options.solver.TimeSolver.movie = 0; % Display movie if 1, last frame otherwise
-options.solver.TimeSolver.save = 0; % 1 --> save during fscan, Warning: memory storage
-
-
-%% Run MOLI
-% Solve the MOLI
-[results,params,options] = MOLI_Control(params,options);
-
-%% Return to HeLaZ workspace
-cd ../../../HeLaZ/wk
diff --git a/matlab/create_gif_1D.m b/matlab/create_gif_1D.m
deleted file mode 100644
index 296117f2103b1dcb7b4c39f5e5f5bc1a86065f39..0000000000000000000000000000000000000000
--- a/matlab/create_gif_1D.m
+++ /dev/null
@@ -1,52 +0,0 @@
-title1 = GIFNAME;
-FIGDIR = BASIC.RESDIR;
-if ~exist(FIGDIR, 'dir')
- mkdir(FIGDIR)
-end
-
-GIFNAME = [FIGDIR, GIFNAME,'.gif'];
-
-% Set colormap boundaries
-hmax = max(max(max(FIELD)));
-hmin = min(min(min(FIELD)));
-flag = 0;
-if hmax == hmin
- disp('Warning : h = hmin = hmax = const')
-else
-% Setup figure frame
-fig = figure;
- plot(X,FIELD(:,1)); % to set up
- axis tight manual % this ensures that getframe() returns a consistent size
- in = 1;
- nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
- for n = FRAMES % loop over selected frames
- scale = max(FIELD(:,n))*SCALING + (1-SCALING);
- plot(X,FIELD(:,n)/scale,linestyle);
- if (YMIN ~= YMAX && XMIN ~= XMAX)
- ylim([YMIN,YMAX]); xlim([XMIN,XMAX]);
- end
- title(['$t \approx$', sprintf('%.3d',ceil(T(n))), ', scaling = ',sprintf('%.1e',scale)]);
- xlabel(XNAME); ylabel(FIELDNAME);
- drawnow
- % Capture the plot as an image
- frame = getframe(fig);
- im = frame2im(frame);
- [imind,cm] = rgb2ind(im,64);
- % Write to the GIF File
- if in == 1
- imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
- else
- imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY);
- end
- % terminal info
- while nbytes > 0
- fprintf('\b')
- nbytes = nbytes - 1;
- end
- nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
- in = in + 1;
- end
- disp(' ')
- disp(['Gif saved @ : ',GIFNAME])
-end
-
diff --git a/matlab/create_gif_1D_phi.m b/matlab/create_gif_1D_phi.m
deleted file mode 100644
index f68761cf852d1cdfcfd2c22edf4b47147235a810..0000000000000000000000000000000000000000
--- a/matlab/create_gif_1D_phi.m
+++ /dev/null
@@ -1,62 +0,0 @@
-title1 = GIFNAME;
-FIGDIR = BASIC.RESDIR;
-if ~exist(FIGDIR, 'dir')
- mkdir(FIGDIR)
-end
-
-GIFNAME = [FIGDIR, GIFNAME,'.gif'];
-
-% Set colormap boundaries
-hmax = max(max(max(FIELD)));
-hmin = min(min(min(FIELD)));
-fieldabsmin = min(min(abs(FIELD(:,:))));
-fieldabsmax = max(max(abs(FIELD(:,:))));
-flag = 0;
-if hmax == hmin
- disp('Warning : h = hmin = hmax = const')
-else
-% Setup figure frame
-fig = figure;
-subplot(211)
- plot(X,FIELD(:,1)); % to set up
- axis tight manual % this ensures that getframe() returns a consistent size
- in = 1;
- nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
- for n = FRAMES % loop over selected frames
- scale = max(FIELD(:,n))*SCALING + (1-SCALING);
-
- subplot(211)
- plot(X,FIELD(:,n)/scale,linestyle);
- if (YMIN ~= YMAX && XMIN ~= XMAX)
- ylim([YMIN,YMAX]); xlim([XMIN,XMAX]);
- end
- title(['$t \approx$', sprintf('%.3d',ceil(T(n))), ', scaling = ',sprintf('%.1e',scale)]);
- xlabel(XNAME); ylabel(FIELDNAME);
- drawnow
- % Capture the plot as an image
- frame = getframe(fig);
- im = frame2im(frame);
- [imind,cm] = rgb2ind(im,64);
- % Write to the GIF File
- if in == 1
- imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
- else
- imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY);
- end
- % terminal info
- while nbytes > 0
- fprintf('\b')
- nbytes = nbytes - 1;
- end
- nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
- in = in + 1;
-
- subplot(212)
- ylim([fieldabsmin,fieldabsmax]); xlim([T(FRAMES(1)),T(FRAMES(end))]); hold on;
- plot(T(n),abs(max(FIELD(:,n))),'.k'); xlabel('$\rho_s/c_s$')
-
- end
- disp(' ')
- disp(['Gif saved @ : ',GIFNAME])
-end
-
diff --git a/matlab/create_gif_3D.m b/matlab/create_gif_3D.m
deleted file mode 100644
index 646080fc0920fddb4f0d48db3d81fa92ccb8cf2d..0000000000000000000000000000000000000000
--- a/matlab/create_gif_3D.m
+++ /dev/null
@@ -1,48 +0,0 @@
-title1 = GIFNAME;
-FIGDIR = BASIC.RESDIR;
-GIFNAME = [FIGDIR, GIFNAME,'.gif'];
-
-% Setup figure frame
-fig = figure('Color','white','Position', [100, 100, 400, 400]);
- plt = @(x) squeeze(reshape(x(:,:,1),[],1));
- scale_x = max(abs(plt(X(:,:,1))));
- scale_y = max(abs(plt(Y(:,:,1))));
- scale_z = max(abs(plt(Z(:,:,1))));
- plot3(plt(X)/scale_x,plt(Y)/scale_y,plt(Z)/scale_z,'.k','MarkerSize',MARKERSIZE);
- view(VIEW);
- axis tight manual % this ensures that getframe() returns a consistent size
- in = 1;
- nbytes = fprintf(2,'frame %d/%d',in,numel(FRAMES));
- for n = FRAMES % loop over selected frames
- plt = @(x) squeeze(reshape(x(:,:,n),[],1));
- scale_x = max(abs(plt(X)));
- scale_y = max(abs(plt(Y)));
- scale_z = max(abs(plt(Z)));
- plot3(plt(X)/scale_x,plt(Y)/scale_y,plt(Z)/scale_z,'.k','MarkerSize',MARKERSIZE);
- view(VIEW);
- xlabel(XNAME); ylabel(YNAME); zlabel(ZNAME);
- title(['$t \approx$', sprintf('%.3d',ceil(T(n)))...
- ,', scaling = ',sprintf('%.1e',scale_x)]);
- grid on;
- drawnow
- % Capture the plot as an image
- frame = getframe(fig);
- im = frame2im(frame);
- [imind,cm] = rgb2ind(im,32);
- % Write to the GIF File
- if in == 1
- imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
- else
- imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY);
- end
- % terminal info
- while nbytes > 0
- fprintf('\b')
- nbytes = nbytes - 1;
- end
- nbytes = fprintf(2,'frame %d/%d',n,numel(FRAMES));
- in = in + 1;
- end
- disp(' ')
- disp(['Gif saved @ : ',GIFNAME])
-
diff --git a/matlab/create_gif_5D.m b/matlab/create_gif_5D.m
deleted file mode 100644
index 3acb6293da6814815df238ec07e268d0ac7c21e8..0000000000000000000000000000000000000000
--- a/matlab/create_gif_5D.m
+++ /dev/null
@@ -1,63 +0,0 @@
-title1 = GIFNAME;
-FIGDIR = BASIC.RESDIR;
-GIFNAME = [FIGDIR, GIFNAME,'.gif'];
-
-sz = size(FIELD);
-
-% Setup figure frame
-fig = figure('Color','white','Position', [100, 100, sz(2)*400, 400]);
- for ij_ = 1:sz(2)
- subplot(100+sz(2)*10+ij_)
-% pclr = imagesc(X,Y,squeeze(FIELD(:,ij_,:,FRAMES(1))));
- pclr = imagesc(X,Y,squeeze(log(FIELD(:,ij_,:,FRAMES(1)))));
- xlabel('$k_r$');
- if ij_ == 1
- ylabel('$P$(max o. $k_z$)');
- else
- yticks([])
- end
- LEGEND = ['$|',FIELDNAME,'^{p',num2str(ij_-1),'}|$']; title(LEGEND);
- end
-% colormap gray
- axis tight manual % this ensures that getframe() returns a consistent size
- in = 1;
- nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,1,:)));
- for n = FRAMES % loop over selected frames
-% scale = max(max(max(abs(FIELD(:,:,:,n)))));
- for ij_ = 1:sz(2)
- subplot(100+sz(2)*10+ij_)
- scale = max(max(max(abs(FIELD(:,ij_,:,n)))));
- pclr = imagesc(X,Y,squeeze(FIELD(:,ij_,:,n))/scale); caxis([0,1]);
-% pclr = imagesc(X,Y,squeeze(log(FIELD(:,ij_,:,n))));
- xlabel(XNAME);
- if ij_ == 1
- ylabel(YNAME);
- else
- yticks([])
- end
- LEGEND = ['$|',FIELDNAME,'^{p',num2str(ij_-1),'}|$'];
- title([LEGEND,', amp = ',sprintf('%.1e',scale)]);
- end
- suptitle(['$t \approx$', sprintf('%.3d',ceil(T(n)))]);
- drawnow
- % Capture the plot as an image
- frame = getframe(fig);
- im = frame2im(frame);
- [imind,cm] = rgb2ind(im,32);
- % Write to the GIF File
- if in == 1
- imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
- else
- imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY);
- end
- % terminal info
- while nbytes > 0
- fprintf('\b')
- nbytes = nbytes - 1;
- end
- nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,1,:)));
- in = in + 1;
- end
- disp(' ')
- disp(['Gif saved @ : ',GIFNAME])
-
diff --git a/matlab/create_gif_imagesc.m b/matlab/create_gif_imagesc.m
deleted file mode 100644
index 78864d803a35a8bba7deb3e08a4e9f70c709a5b5..0000000000000000000000000000000000000000
--- a/matlab/create_gif_imagesc.m
+++ /dev/null
@@ -1,49 +0,0 @@
-title1 = GIFNAME;
-FIGDIR = BASIC.RESDIR;
-GIFNAME = [FIGDIR, GIFNAME,'.gif'];
-
-% Set colormap boundaries
-hmax = max(max(max(FIELD)));
-hmin = min(min(min(FIELD)));
-
-flag = 0;
-if hmax == hmin
- disp('Warning : h = hmin = hmax = const')
-else
-% Setup figure frame
-fig = figure('Color','white','Position', [100, 100, 400, 400]);
- imagesc(X,Y,FIELD(:,:,1)); % to set up
- colormap gray
- axis tight manual % this ensures that getframe() returns a consistent size
- in = 1;
- nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
- for n = FRAMES % loop over selected frames
- scale = max(max(abs(FIELD(:,:,n))));
- pclr = imagesc(X,Y,FIELD(:,:,n)/scale); % frame plot
-% caxis([min(min(FIELD(:,:,n))),max(max(FIELD(:,:,n)))]);
- xlabel(XNAME); ylabel(YNAME); %colorbar;
- title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
- ,', scaling = ',sprintf('%.1e',scale)]);
- drawnow
- % Capture the plot as an image
- frame = getframe(fig);
- im = frame2im(frame);
- [imind,cm] = rgb2ind(im,32);
- % Write to the GIF File
- if in == 1
- imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
- else
- imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY);
- end
- % terminal info
- while nbytes > 0
- fprintf('\b')
- nbytes = nbytes - 1;
- end
- nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
- in = in + 1;
- end
- disp(' ')
- disp(['Gif saved @ : ',GIFNAME])
-end
-
diff --git a/matlab/create_gif_surf.m b/matlab/create_gif_surf.m
deleted file mode 100644
index 087b788c87053b3bb3822e41e27b84e4521f3e94..0000000000000000000000000000000000000000
--- a/matlab/create_gif_surf.m
+++ /dev/null
@@ -1,56 +0,0 @@
-title1 = GIFNAME;
-FIGDIR = BASIC.RESDIR;
-GIFNAME = [FIGDIR, GIFNAME,'.gif'];
-
-% Set colormap boundaries
-hmax = max(max(max(FIELD)));
-hmin = min(min(min(FIELD)));
-
-flag = 0;
-if hmax == hmin
- disp('Warning : h = hmin = hmax = const')
-else
-% Setup figure frame
-fig = figure('Color','white','Position', [100, 100, 400, 400]);
- surf(X,Y,FIELD(:,:,1)); % to set up
- colormap gray
- axis tight manual % this ensures that getframe() returns a consistent size
- if INTERP
- shading interp;
- end
- in = 1;
- nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
- for n = FRAMES % loop over selected frames
- scale = max(max(abs(FIELD(:,:,n))));
- pclr = surf(X,Y,FIELD(:,:,n)/scale); % frame plot
- if INTERP
- shading interp;
- end
- set(pclr, 'edgecolor','none'); axis square;
-% caxis([min(min(FIELD(:,:,n))),max(max(FIELD(:,:,n)))]);
- xlabel(XNAME); ylabel(YNAME); %colorbar;
- title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
- ,', scaling = ',sprintf('%.1e',scale)]);
- drawnow
- % Capture the plot as an image
- frame = getframe(fig);
- im = frame2im(frame);
- [imind,cm] = rgb2ind(im,32);
- % Write to the GIF File
- if in == 1
- imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
- else
- imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY);
- end
- % terminal info
- while nbytes > 0
- fprintf('\b')
- nbytes = nbytes - 1;
- end
- nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
- in = in + 1;
- end
- disp(' ')
- disp(['Gif saved @ : ',GIFNAME])
-end
-
diff --git a/matlab/half_2_full_cc_2D.m b/matlab/half_2_full_cc_2D.m
deleted file mode 100644
index 481ee1620447f1e502814b78366b27213284aece..0000000000000000000000000000000000000000
--- a/matlab/half_2_full_cc_2D.m
+++ /dev/null
@@ -1,27 +0,0 @@
-function [ FF_f ] = half_2_full_cc_2D( FF_h )
-%half_2_full_cc_2D Retrieve full domain frequency from half domain
-% Meant for complex conjugate symmetric fields
-
-[Nkr,Nkz] = size(FF_h);
-
-if Nkr > Nkz
- FF_f = zeros(max(size(FF_h)));
- FF_f(1:Nkr,1:Nkz) = FF_h;
-
- for ikr = 1:Nkr
- for ikz = Nkr/2+2:Nkr
- FF_f(ikr,ikz) = FF_f(Nkr-ikr+1,Nkr-ikz+1);
- end
- end
-else
- FF_f = zeros(max(size(FF_h)));
- FF_f(1:Nkr,1:Nkz) = FF_h;
-
- for ikz = 1:Nkz
- for ikr = Nkz/2+2:Nkz
- FF_f(ikr,ikz) = FF_f(Nkz-ikr+1,Nkz-ikz+1);
- end
- end
-end
-end
-
diff --git a/matlab/helaz_analysis.m b/matlab/helaz_analysis.m
deleted file mode 100644
index e623ae83b73cfaf865ffb667d77c99ea8eeb00c8..0000000000000000000000000000000000000000
--- a/matlab/helaz_analysis.m
+++ /dev/null
@@ -1,27 +0,0 @@
-%% HeLaZ data
-filename = 'results_00.h5';
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Load the data
-moment = 'Ni00';
-
-kr = h5read(filename,['/data/var2d/' moment '/coordkr']);
-kz = h5read(filename,['/data/var2d/' moment '/coordkz']);
-timeNi = h5read(filename,'/data/var2d/time');
-Nipj = zeros(numel(timeNi),numel(kr),numel(kz));
-for it = 1:numel(timeNi)
- tmp = h5read(filename,['/data/var2d/', moment,'/', num2str(it,'%06d')]);
- Nipj(it,:,:) = tmp.real + 1i * tmp.imaginary;
-end
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Plot growth rate vs kz
-K_RICCI = 1; %% add a sqrt(1+tau) to the kperps
-ikr = 1; %% Fix the kr value
-plot_gamma_vs_k;
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Plot Ni00 evolution
-ikr = 1; %% Fix the kr value
-plot_Ni00_t_evolution;
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
diff --git a/matlab/plot_Ni00_t_evolution.m b/matlab/plot_Ni00_t_evolution.m
deleted file mode 100644
index d5aa35dcc29fb37919669b6265675e2f8b1eaa7d..0000000000000000000000000000000000000000
--- a/matlab/plot_Ni00_t_evolution.m
+++ /dev/null
@@ -1,44 +0,0 @@
-%% Plot the time evolution of the firt ion moment
-default_plots_options % Script to set up default plot variables
-
-fig = figure;
-
-LEGEND = [];
-
-x1 = timeNi;
-
-for ikz = 1:2:numel(kz)
-
- linename = ['$k_r = ',num2str(kr(ikr)),'$, ','$k_z = ',num2str(kz(ikz)),'$'];
- y1 = abs(Nipj(:,ikr,ikz));
- semilogy(x1,y1,'DisplayName',linename)
- LEGEND = [LEGEND, linename];
- hold on
-
-end
-
-for ikz = 1:2:numel(kz)
-
- semilogy(x1(itmin:end),...
- exp(gammas(ikr,ikz)*x1(itmin:end) + shifts(ikr,ikz)),...
- 'Color', 'k', 'LineStyle', '--','HandleVisibility','off')
-
-end
-
-LEGEND = [LEGEND, 'fits'];
-
-TITLE = [];
-TITLE = [TITLE,'$\eta_n=',num2str(MODEL.eta_n),'$, '];
-TITLE = [TITLE,'$\eta_B=',num2str(MODEL.eta_B),'$, '];
-TITLE = [TITLE, '$\nu=',num2str(MODEL.nu),'$, '];
-TITLE = [TITLE, '$(P,J)=(',num2str(GRID.pmaxe),',',num2str(GRID.jmaxe),')$'];
-title(TITLE);
-grid on
-xlabel('$t$')
-ylabel(['$|',moment,'|$'])
-
-%% Saving fig
-FIGNAME = 'Ni00_t';
-if SAVEFIG
- save_figure;
-end
\ No newline at end of file
diff --git a/matlab/plot_gamma_vs_k.m b/matlab/plot_gamma_vs_k.m
deleted file mode 100644
index 250b4a1246d4b1dc5974cbb9fccaf44becc459e8..0000000000000000000000000000000000000000
--- a/matlab/plot_gamma_vs_k.m
+++ /dev/null
@@ -1,45 +0,0 @@
-%Plot growth rate vs kz
-default_plots_options % Script to set up default plot variables
-% with a linear fit of the log evolution
-gammas = zeros(numel(kr),numel(kz));
-shifts = zeros(numel(kr),numel(kz));
-
-if K_RICCI
- factor = sqrt(1+MODEL.tau_i);
- fchar = '\times(1+\tau)^{1/2}$';
-else
- factor = 1;
- fchar = '$';
-end
-
-% Linear fit of log(Napj)
-x1 = timeNi;
-itmin = ceil(0.9 * numel(timeNi)); %Take a subset of the time evolution
-
-for ikz = 1:numel(kz)
- fit = polyfit(x1(itmin:end),log(abs(Nipj(itmin:end,ikr,ikz))),1);
- gammas(ikr,ikz) = fit(1);
- shifts(ikr,ikz) = fit(2);
-end
-
-fig = figure;
-linename = ['$(P,J)=(',num2str(GRID.pmaxe),',',num2str(GRID.jmaxe),')$'];
-plot(factor*kz,gammas(ikr,:),'DisplayName',linename);
-
-TITLE = [];
-TITLE = [TITLE,'$\eta_n=',num2str(1.0/MODEL.eta_n),'$, '];
-TITLE = [TITLE,'$\eta_B=',num2str(MODEL.eta_B),'$, '];
-TITLE = [TITLE, '$\nu=',num2str(MODEL.nu),'$, '];
-%TITLE = [TITLE, '$k_z=',num2str(GRID.kz),'$'];
-
-title(TITLE);
-grid on
-legend('show')
-xlabel(['$k_z',fchar])
-ylabel('$\gamma L_\perp/c_{s} $')
-
-%% Saving fig
-if SAVEFIG
- FIGNAME = 'gamma_k';
- save_figure;
-end
diff --git a/wk/analysis_1D.m b/wk/analysis_1D.m
deleted file mode 100644
index bcfffa5749e7677f70cfb5ea251b861a6b38f6fe..0000000000000000000000000000000000000000
--- a/wk/analysis_1D.m
+++ /dev/null
@@ -1,233 +0,0 @@
-default_plots_options
-%% load resulTs2D
-JOBNUM = 00;
-load_results
-Ni00 = squeeze(Ni00);
-Ne00 = squeeze(Ne00);
-PHI = squeeze(PHI);
-Ts2D = Ts2D';
-Ns = numel(Ts2D);
-dt_samp = mean(diff(Ts2D));
-% Build grids
-Nkx = numel(kx); Nky = numel(ky);
-[ky,kx] = meshgrid(ky,kx);
-Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky);
-dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1);
-Lk = max(Lkx,Lky);
-
-dr = 2*pi/Lk; dz = 2*pi/Lk;
-Nx = max(Nkx,Nky) * (Nkx > 1) + (Nkx == 1); Ny = Nky;
-r = dr*(-Nx/2:(Nx/2-1)); Lx = max(r)-min(r);
-z = dz*(-Ny/2:(Ny/2-1)); Ly = max(z)-min(z);
-[YY,XX] = meshgrid(z,r);
-% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% IFFT
-ne00 = zeros(Ny, Ns);
-ni00 = zeros(Ny, Ns);
-phi = zeros(Ny, Ns);
-
-for it = 1:numel(Ts2D)
- NE_ = Ne00(:,it); NN_ = Ni00(:,it); PH_ = PHI(:,it);
- F_ = (ifft((NE_),Ny));
- ne00(:,it)= real(fftshift(F_));
- F_ = (ifft((NN_),Ny));
- ni00(:,it) = real(fftshift(F_));
- F_ = (ifft((PH_),Ny));
- phi(:,it) = real(fftshift(F_));
-end
-
-% Post processing
-ne_00 = zeros(1,Ns); % Time evolution of ne(r,z) at origin
-ni_00 = zeros(1,Ns); % .
-phi_00 = zeros(1,Ns); % .
-E_pot = zeros(1,Ns); % Potential energy n^2
-E_kin = zeros(1,Ns); % Kinetic energy grad(phi)^2
-ExB = zeros(1,Ns); % ExB drift intensity \propto |\grad \phi|
-CFL = zeros(1,Ns); % CFL time step
-Ddz = 1i*ky;
-[~,iz0] = min(abs(z)); % index of z==0
-[~,iky1] = min(abs(ky-round(1/dky)*dky)); % index of ky==1
-for it = 1:numel(Ts2D)
- NE_ = squeeze(Ne00(:,it)); NN_ = squeeze(Ni00(:,it)); PH_ = squeeze(PHI(:,it));
-
- ne_00(it) = ne00(iz0,it);
-
- ni_00(it) = ni00(iz0,it);
-
- phi_00(it) = phi(iz0,it);
-
- E_pot(it) = pi/Ly*sum(abs(NN_).^2)/Nky; % integrate through Parseval id
- E_kin(it) = pi/Ly*sum(abs(Ddz.*PH_).^2)/Nky;
- ExB(it) = max(abs(phi(3:end,it)-phi(1:end-2,it))'/(2*dz));
-end
-E_kin_KY = abs(Ddz.*PHI(:,it)).^2;
-dEdt = diff(E_pot+E_kin)./diff(Ts2D);
-
-%% Compute growth rate
-gamma_Ne = zeros(1,Nky);
-gamma_Ni = zeros(1,Nky);
-gamma_PH = zeros(1,Nky);
-tend = Ts2D(end); tstart = 0.6*tend;
-plt = @(x) squeeze(abs(x));
-for iky = 1:Nky
- gamma_Ne(iky) = LinearFit_s(Ts2D,plt(Ne00(iky,:)),tstart,tend);
- gamma_Ni(iky) = LinearFit_s(Ts2D,plt(Ni00(iky,:)),tstart,tend);
- gamma_PH(iky) = LinearFit_s(Ts2D,plt(PHI(iky,:)),tstart,tend);
-end
-gamma_Ne = real(gamma_Ne .* (gamma_Ne>=0.0));
-gamma_Ni = real(gamma_Ni .* (gamma_Ni>=0.0));
-gamma_PH = real(gamma_PH .* (gamma_PH>=0.0));
-
-%% PLOTs2D
-if 0
-%% Time evolutions
-fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM)];
- subplot(221)
- semilogy(Ts2D,abs(ne_00),'-','DisplayName','$n_e^{00}$'); hold on;
- semilogy(Ts2D,abs(ni_00),'-','DisplayName','$n_i^{00}$');
- grid on; xlabel('$t$'); ylabel('$|n_a(x=0,y=0)|$');
- subplot(222)
- semilogy(Ts2D,abs(Ni_gm),'-','DisplayName','$\phi$')
- grid on; xlabel('$t$'); ylabel('$|\tilde n(k_r\approx 1,k_z\approx 1)|$');
- subplot(223)
- semilogy(Ts2D,E_kin+E_pot,'-','DisplayName','$\sum|ik\tilde\phi_i|^2+\sum|\tilde n_i|^2$')
- hold on;
- grid on; xlabel('$t$'); ylabel('$E$'); legend('show');
-FMT = '.fig'; save_figure
-end
-if 1
-%% Growth rate
-fig = figure; FIGNAME = ['growth_rate',sprintf('_%.2d',JOBNUM)];
-plt = @(x) circshift(x,Nky/2-1);
- subplot(221)
- plot(plt(ky),plt(gamma_Ne),'-'); hold on; xlim([0.0,max(ky)+dky]);
- grid on; xlabel('$k_z$'); ylabel('$\gamma(N_e^{00})$');
- subplot(222)
- plot(plt(ky),plt(gamma_Ni),'-'); hold on; xlim([0.0,max(ky)+dky]);
- grid on; xlabel('$k_z$'); ylabel('$\gamma(N_i^{00})$');
- subplot(223)
- plot(plt(ky),plt(gamma_PH),'-'); hold on; xlim([0.0,max(ky)+dky]);
- grid on; xlabel('$k_z$'); ylabel('$\gamma(\tilde\phi)$'); legend('show');
-FMT = '.fig'; save_figure
-end
-%% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-t0 = 0;
-skip_ = 1;
-DELAY = 0.01*skip_;
-FRAMES = floor(t0/dt_samp)+1:skip_:numel(Ts2D);
-linestyle = 'o-.';
-if 0
-%% Density electron
-GIFNAME = ['ne',sprintf('_%.2d',JOBNUM)];
-FIELD = real(ne00); X = z; T = Ts2D;
-FIELDNAME = '$n_e^{00}/\max(n_e^{00})$'; XNAME = '$z$';
-XMIN = -L/2-2; XMAX = L/2+1; YMIN = -1.1; YMAX = 1.1;
-create_gif_1D
-%% Density ion
-GIFNAME = ['ni',sprintf('_%.2d',JOBNUM)];
-FIELD = real(ni00); X = z; T = Ts2D;
-FIELDNAME = '$n_i^{00}/\max(n_i^{00})$'; XNAME = '$z$';
-XMIN = -L/2-2; XMAX = L/2+1; YMIN = -1.1; YMAX = 1.1;
-create_gif_1D
-%% Phi
-GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)];
-FIELD = real(phi); X = z; T = Ts2D;
-FIELDNAME = '$\phi/\max(\phi)$'; XNAME = '$z$';
-XMIN = -L/2-2; XMAX = L/2+1; YMIN = -1.1; YMAX = 1.1;
-create_gif_1D
-%% Density electron frequency
-GIFNAME = ['Ni',sprintf('_%.2d',JOBNUM)];
-FIELD = (abs(Ni00)); X = (ky); T = Ts2D;
-FIELDNAME = '$|N_i^{00}|$'; XNAME = '$k_z$';
-XMIN = min(ky)-0.5; XMAX = max(ky)+.5; YMIN = -0.1; YMAX = 1.1;
-create_gif_1D
-end
-
-if 0
-%% Space-Time diagrams at r = 0
-plt = @(x) real(x);
-fig = figure; FIGNAME = ['z_space_time_diag',sprintf('_%.2d',JOBNUM)];
- [TY,TX] = meshgrid(Ts2D,z);
-subplot(221)% density
- pclr = pcolor(TX,TY,(plt(ne00))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$z\,(r=0)$'); ylabel('$t$'); title('$n_e^{00}$');
-subplot(222)% density
- pclr = pcolor(TX,TY,(plt(ni00))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$z\,(r=0)$'); ylabel('$t$'); title('$n_i^{00}$');
-subplot(223)% density
- pclr = pcolor(TX,TY,(plt(phi))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$z\,(r=0)$'); ylabel('$t$'); title('$\phi$');
-FMT = '.fig'; save_figure
-
-%% Space-Time diagrams at kx=0
-plt = @(x) abs(x);
-fig = figure; FIGNAME = ['ky_space_time_diag',sprintf('_%.2d',JOBNUM)];
- [TY,TX] = meshgrid(Ts2D,ky);
-subplot(221)% density
- pclr = pcolor(TX,TY,(plt(Ne00))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$ky$'); ylabel('$t$'); title('$\textrm{Re}(N_e^{00})|_{k_r=0}$');
-subplot(222)% density
- pclr = pcolor(TX,TY,(plt(Ni00))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$ky$'); ylabel('$t$'); title('$\textrm{Re}(N_i^{00})|_{k_r=0}$');
-subplot(223)% density
- pclr = pcolor(TX,TY,(plt(PH))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$ky$'); ylabel('$t$'); title('$\textrm{Re}(\tilde\phi)|_{k_r=0}$');
-FMT = '.fig'; save_figure
-end
-
-if 0
-%% Mode time evolution
-[~,ik05] = min(abs(ky-0.50));
-[~,ik10] = min(abs(ky-1.00));
-[~,ik15] = min(abs(ky-1.50));
-[~,ik20] = min(abs(ky-2.00));
-
-fig = figure; FIGNAME = ['frame',sprintf('_%.2d',JOBNUM)];
- subplot(221); plt = @(x) abs(squeeze(x));
- semilogy(Ts2D,plt(PHI(ik05,:)),'DisplayName', '$k_z = 0.5$'); hold on
- semilogy(Ts2D,plt(PHI(ik10,:)),'DisplayName', '$k_z = 1.0$')
- semilogy(Ts2D,plt(PHI(ik15,:)),'DisplayName', '$k_z = 1.5$')
- semilogy(Ts2D,plt(PHI(ik20,:)),'DisplayName', '$k_z = 2.0$')
- xlabel('$t$'); ylabel('$|\tilde\phi|$'); legend('show');
- semilogy(Ts2D,plt(PHI(ik05,end)).*exp(gamma_PH(ik05).*(Ts2D-Ts2D(end))),'--k')
- semilogy(Ts2D,plt(PHI(ik10,end)).*exp(gamma_PH(ik10).*(Ts2D-Ts2D(end))),'--k')
- plot(tstart*[1 1],ylim,'k-','LineWidth',0.5);
- plot(tend*[1 1],ylim,'k-','LineWidth',0.5);
- subplot(222); plt = @(x) squeeze(abs(x));
- semilogy(Ts2D,plt(Ni00(ik05,:)),'DisplayName', '$k_z = 0.5$'); hold on
- semilogy(Ts2D,plt(Ni00(ik10,:)),'DisplayName', '$k_z = 1.0$')
- semilogy(Ts2D,plt(Ni00(ik15,:)),'DisplayName', '$k_z = 1.5$')
- semilogy(Ts2D,plt(Ni00(ik20,:)),'DisplayName', '$k_z = 2.0$')
- xlabel('$t$'); ylabel('$|N_i^{00}|$'); legend('show');
- semilogy(Ts2D,plt(Ni00(ik05,end)).*exp(gamma_Ni(ik05).*(Ts2D-Ts2D(end))),'--k')
- semilogy(Ts2D,plt(Ni00(ik10,end)).*exp(gamma_Ni(ik10).*(Ts2D-Ts2D(end))),'--k')
- plot(tstart*[1 1],ylim,'k-','LineWidth',0.5);
- plot(tend*[1 1],ylim,'k-','LineWidth',0.5);
- subplot(223); plt = @(x) squeeze(abs(x));
- semilogy(Ts2D,plt(Ne00(ik05,:)),'DisplayName', '$k_z = 0.5$'); hold on
- semilogy(Ts2D,plt(Ne00(ik10,:)),'DisplayName', '$k_z = 1.0$')
- semilogy(Ts2D,plt(Ne00(ik15,:)),'DisplayName', '$k_z = 1.5$')
- semilogy(Ts2D,plt(Ne00(ik20,:)),'DisplayName', '$k_z = 2.0$')
- xlabel('$t$'); ylabel('$|N_e^{00}|$'); legend('show');
- semilogy(Ts2D,plt(Ne00(ik05,end)).*exp(gamma_Ne(ik05).*(Ts2D-Ts2D(end))),'--k')
- semilogy(Ts2D,plt(Ne00(ik10,end)).*exp(gamma_Ne(ik10).*(Ts2D-Ts2D(end))),'--k')
- plot(tstart*[1 1],ylim,'k-','LineWidth',0.5);
- plot(tend*[1 1],ylim,'k-','LineWidth',0.5);
-FMT = '.fig'; save_figure
-end
-
-if 0
-%% Show frame
-it = min(70,numel(Ts2D));
-fig = figure; FIGNAME = ['frame',sprintf('_%.2d',JOBNUM)];
- subplot(221); plt = @(x) (abs(x));
- plot(ky,plt(PH(:,it)))
- xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$\hat\phi$');
- subplot(222); plt = @(x) (abs(x));
- plot(ky,plt(Ni00(:,it)))
- xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$\hat n_i^{00}$');
- subplot(223); plt = @(x) (abs(x));
- plot(ky,plt(Ne00(:,it)))
- xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$\hat n_e^{00}$');
-FMT = '.fig'; save_figure
-end
\ No newline at end of file
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
deleted file mode 100644
index 50a40384355df85c07257905a776b7e2e718ef42..0000000000000000000000000000000000000000
--- a/wk/analysis_2D.m
+++ /dev/null
@@ -1,599 +0,0 @@
-addpath(genpath('../matlab')) % ... add
-for i_ = 1
-% for ETA_ =[0.6:0.1:0.9]
-%% Load results
-if 0% Local results
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='v2.8_kobayashi/100x50_L_50_P_6_J_3_eta_0.71429_nu_1e-02_PAGK_CLOS_0_mu_0e+00';
- BASIC.RESDIR = ['../results/',outfile,'/'];
- BASIC.MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
- CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
- system(CMD);
-end
-if 1% Marconi results
-outfile ='';
-outfile ='';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/out.txt';
-% outfile = outcl{i_};
-% load_marconi(outfile);
-BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
-BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
-end
-
-%%
-% JOBNUM = 1; load_results;
-JOBNUMMAX = 20; compile_results %Compile the results from first output found to JOBNUMMAX if existing
-
-%% Retrieving max polynomial degree and sampling info
-Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe);
-Npi = numel(Pi); Nji = numel(Ji); [JI,PI] = meshgrid(Ji,Pi);
-Ns5D = numel(Ts5D);
-Ns2D = numel(Ts2D);
-% renaming and reshaping quantity of interest
-Ts5D = Ts5D';
-Ts2D = Ts2D';
-
-%% Build grids
-Nkx = numel(kx); Nky = numel(ky);
-[ky,kx] = meshgrid(ky,kx);
-Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky);
-dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1);
-KPERP2 = ky.^2+kx.^2;
-[~,ikx0] = min(abs(kx)); [~,iky0] = min(abs(ky));
-
-Lk = max(Lkx,Lky);
-dr = 2*pi/Lk; dz = 2*pi/Lk;
-Nx = max(Nkx,Nky); Ny = Nx;
-r = dr*(-Nx/2:(Nx/2-1)); Lx = max(r)-min(r);
-z = dz*(-Ny/2:(Ny/2-1)); Ly = max(z)-min(z);
-[ZZ,RR] = meshgrid(z,r);
-
-%% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-disp('Analysis :')
-disp('- iFFT')
-% IFFT (Lower case = real space, upper case = frequency space)
-ne00 = zeros(Nx,Ny,Ns2D); % Gyrocenter density
-ni00 = zeros(Nx,Ny,Ns2D);
-dzTe = zeros(Nx,Ny,Ns2D);
-dzTi = zeros(Nx,Ny,Ns2D);
-dzni = zeros(Nx,Ny,Ns2D);
-np_i = zeros(Nx,Ny,Ns5D); % Ion particle density
-si00 = zeros(Nx,Ny,Ns5D);
-phi = zeros(Nx,Ny,Ns2D);
-dens_e = zeros(Nx,Ny,Ns2D);
-dens_i = zeros(Nx,Ny,Ns2D);
-temp_e = zeros(Nx,Ny,Ns2D);
-temp_i = zeros(Nx,Ny,Ns2D);
-drphi = zeros(Nx,Ny,Ns2D);
-dzphi = zeros(Nx,Ny,Ns2D);
-dr2phi = zeros(Nx,Ny,Ns2D);
-
-for it = 1:numel(Ts2D)
- NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it);
- ne00(:,:,it) = real(fftshift(ifft2((NE_),Nx,Ny)));
- ni00(:,:,it) = real(fftshift(ifft2((NI_),Nx,Ny)));
- phi (:,:,it) = real(fftshift(ifft2((PH_),Nx,Ny)));
- drphi(:,:,it) = real(fftshift(ifft2(1i*kx.*(PH_),Nx,Ny)));
- dr2phi(:,:,it)= real(fftshift(ifft2(-kx.^2.*(PH_),Nx,Ny)));
- dzphi(:,:,it) = real(fftshift(ifft2(1i*ky.*(PH_),Nx,Ny)));
- if(W_DENS && W_TEMP)
- DENS_E_ = DENS_E(:,:,it); DENS_I_ = DENS_I(:,:,it);
- TEMP_E_ = TEMP_E(:,:,it); TEMP_I_ = TEMP_I(:,:,it);
- dzni(:,:,it) = real(fftshift(ifft2(1i*ky.*(DENS_I_),Nx,Ny)));
- dzTe(:,:,it) = real(fftshift(ifft2(1i*ky.*(TEMP_E_),Nx,Ny)));
- dzTi(:,:,it) = real(fftshift(ifft2(1i*ky.*(TEMP_I_),Nx,Ny)));
- dens_e (:,:,it) = real(fftshift(ifft2((DENS_E_),Nx,Ny)));
- dens_i (:,:,it) = real(fftshift(ifft2((DENS_I_),Nx,Ny)));
- temp_e (:,:,it) = real(fftshift(ifft2((TEMP_E_),Nx,Ny)));
- temp_i (:,:,it) = real(fftshift(ifft2((TEMP_I_),Nx,Ny)));
- end
-end
-
-% Building a version of phi only 5D sampling times
-PHI_Ts5D = zeros(Nkx,Nky,Ns5D);
-err = 0;
-for it = 1:numel(Ts5D) % Loop over 5D arrays
- [shift, it2D] = min(abs(Ts2D-Ts5D(it)));
- if shift > err; err = shift; end;
- PHI_Ts5D(:,:,it) = PHI(:,:,it2D);
-end
-if err > 0; disp('WARNING Ts2D and Ts5D are shifted'); end;
-
-Np_i = zeros(Nkx,Nky,Ns5D); % Ion particle density in Fourier space
-
-for it = 1:numel(Ts5D)
- [~, it2D] = min(abs(Ts2D-Ts5D(it)));
- Np_i(:,:,it) = 0;
- for ij = 1:Nji
- Kn = (KPERP2/2.).^(ij-1) .* exp(-KPERP2/2)/(factorial(ij-1));
- Np_i(:,:,it) = Np_i(:,:,it) + Kn.*squeeze(Nipj(1,ij,:,:,it));
- end
-
- np_i(:,:,it) = real(fftshift(ifft2(squeeze(Np_i(:,:,it)),Nx,Ny)));
-end
-
-% Post processing
-disp('- post processing')
-% gyrocenter and particle flux from real space
-GFlux_ri = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ni drphi>
-GFlux_zi = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ni dzphi>
-GFlux_re = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ne drphi>
-GFlux_ze = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ne dzphi>
-PFlux_ri = zeros(1,Ns5D); % Particle flux
-% gyrocenter and particle flux from fourier coefficients
-GFLUX_RI = real(squeeze(sum(sum(-1i*ky.*Ni00.*conj(PHI),1),2)))*(2*pi/Nx/Ny)^2;
-PFLUX_RI = real(squeeze(sum(sum(-1i*ky.*Np_i.*conj(PHI_Ts5D),1),2)))*(2*pi/Nx/Ny)^2;
-% Heat flux
-Q_RI = -squeeze(mean(mean(dzphi.*temp_i,1),2))';
-% Hermite energy spectrum
-epsilon_e_pj = zeros(Pe_max,Je_max,Ns5D);
-epsilon_i_pj = zeros(Pi_max,Ji_max,Ns5D);
-
-phi_maxr_maxz = zeros(1,Ns2D); % Time evol. of the norm of phi
-phi_avgr_maxz = zeros(1,Ns2D); % Time evol. of the norm of phi
-phi_maxr_avgz = zeros(1,Ns2D); % Time evol. of the norm of phi
-phi_avgr_avgz = zeros(1,Ns2D); % Time evol. of the norm of phi
-
-shear_maxr_maxz = zeros(1,Ns2D); % Time evol. of the norm of shear
-shear_avgr_maxz = zeros(1,Ns2D); % Time evol. of the norm of shear
-shear_maxr_avgz = zeros(1,Ns2D); % Time evol. of the norm of shear
-shear_avgr_avgz = zeros(1,Ns2D); % Time evol. of the norm of shear
-
-Ne_norm = zeros(Pe_max,Je_max,Ns5D); % Time evol. of the norm of Napj
-Ni_norm = zeros(Pi_max,Ji_max,Ns5D); % .
-
-Ddr = 1i*kx; Ddz = 1i*ky; lapl = Ddr.^2 + Ddz.^2;
-% Kperp spectrum interpolation
-%full kperp points
-kperp = reshape(sqrt(kx.^2+ky.^2),[numel(kx),1]);
-% interpolated kperps
-nk_noAA = floor(2/3*numel(kx));
-kp_ip = kx;
-[thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip);
-[xn,yn] = pol2cart(thg,rg);
-[ky_s, sortIdx] = sort(ky);
-[xc,yc] = meshgrid(ky_s,kx);
-phi_kp_t = zeros(numel(kp_ip),Ns2D);
-%
-for it = 1:numel(Ts2D) % Loop over 2D arrays
- NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it);
- phi_maxr_maxz(it) = max( max(squeeze(phi(:,:,it))));
- phi_avgr_maxz(it) = max(mean(squeeze(phi(:,:,it))));
- phi_maxr_avgz(it) = mean( max(squeeze(phi(:,:,it))));
- phi_avgr_avgz(it) = mean(mean(squeeze(phi(:,:,it))));
-
- shear_maxr_maxz(it) = max( max(squeeze(-(dr2phi(:,:,it)))));
- shear_avgr_maxz(it) = max(mean(squeeze(-(dr2phi(:,:,it)))));
- shear_maxr_avgz(it) = mean( max(squeeze(-(dr2phi(:,:,it)))));
- shear_avgr_avgz(it) = mean(mean(squeeze(-(dr2phi(:,:,it)))));
-
- GFlux_ri(it) = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lx/Ly;
- GFlux_zi(it) = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lx/Ly;
- GFlux_re(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lx/Ly;
- GFlux_ze(it) = sum(sum(-ne00(:,:,it).*drphi(:,:,it)))*dr*dz/Lx/Ly;
-
- Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,it))).^2,3)),xn,yn);
- phi_kp_t(:,it) = mean(Z_rth,2);
-end
-%
-for it = 1:numel(Ts5D) % Loop over 5D arrays
- [~, it2D] = min(abs(Ts2D-Ts5D(it)));
- Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkx/Nky;
- Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkx/Nky;
- epsilon_e_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nepj(:,:,:,:,it)).^2,3),4);
- epsilon_i_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nipj(:,:,:,:,it)).^2,3),4);
- % Particle flux
- PFlux_ri(it) = sum(sum(np_i(:,:,it).*dzphi(:,:,it2D)))*dr*dz/Lx/Ly;
-end
-
-%% Compute primary instability growth rate
-disp('- growth rate')
-% Find max value of transport (end of linear mode)
-[tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2);
-[~,itmax] = min(abs(Ts2D-tmax));
-tstart = 0.1 * Ts2D(itmax); tend = 0.5 * Ts2D(itmax);
-[~,its2D_lin] = min(abs(Ts2D-tstart));
-[~,ite2D_lin] = min(abs(Ts2D-tend));
-
-g_I = zeros(Nkx,Nky);
-for ikx = 1:Nkx
- for iky = 1:Nky
- [g_I(ikx,iky), ~] = LinearFit_s(Ts2D(its2D_lin:ite2D_lin),squeeze(abs(Ni00(ikx,iky,its2D_lin:ite2D_lin))));
- end
-end
-[gmax_I,ikmax_I] = max(g_I(1,:));
-kmax_I = abs(ky(ikmax_I));
-Bohm_transport = ETAB/ETAN*gmax_I/kmax_I^2;
-
-%% Compute secondary instability growth rate
-disp('- growth rate')
-% Find max value of transport (end of linear mode)
-% [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2);
-% [~,itmax] = min(abs(Ts2D-tmax));
-% tstart = Ts2D(itmax); tend = 1.5*Ts2D(itmax);
-[~,its2D_lin] = min(abs(Ts2D-tstart));
-[~,ite2D_lin] = min(abs(Ts2D-tend));
-
-g_II = zeros(Nkx,Nky);
-for ikx = 1:Nkx
- for iky = 1
- [g_II(ikx,iky), ~] = LinearFit_s(Ts2D(its2D_lin:ite2D_lin),squeeze(abs(Ni00(ikx,iky,its2D_lin:ite2D_lin))));
- end
-end
-[gmax_II,ikmax_II] = max(g_II(1,:));
-kmax_II = abs(kx(ikmax_II));
-
-%% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-default_plots_options
-disp('Plots')
-FMT = '.fig';
-
-if 1
-%% Time evolutions and growth rate
-fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM),'_',PARAMS];
-set(gcf, 'Position', [100, 100, 900, 800])
-subplot(111);
- suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
- ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
- ' $\mu_{hd}=$',num2str(MU)]);
- subplot(421);
- for ip = 1:Pe_max
- for ij = 1:Je_max
- plt = @(x) squeeze(x(ip,ij,:));
- plotname = ['$N_e^{',num2str(ip-1),num2str(ij-1),'}$'];
- clr = line_colors(min(ip,numel(line_colors(:,1))),:);
- lstyle = line_styles(min(ij,numel(line_styles)));
- semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname,...
- 'Color',clr,'LineStyle',lstyle{1}); hold on;
- end
- end
- grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$');
- subplot(423)
- for ip = 1:Pi_max
- for ij = 1:Ji_max
- plt = @(x) squeeze(x(ip,ij,:));
- plotname = ['$N_i^{',num2str(ip-1),num2str(ij-1),'}$'];
- clr = line_colors(min(ip,numel(line_colors(:,1))),:);
- lstyle = line_styles(min(ij,numel(line_styles)));
- semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname,...
- 'Color',clr,'LineStyle',lstyle{1}); hold on;
- end
- end
- grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); xlabel('$t c_s/R$')
- subplot(222)
- plot(Ts0D,GGAMMA_RI*(2*pi/Nx/Ny)^2); hold on;
-% plot(Ts2D,GFLUX_RI)
- plot(Ts0D,PGAMMA_RI*(2*pi/Nx/Ny)^2);
-% plot(Ts5D,PFLUX_RI,'--');
- legend(['Gyro. flux';'Part. flux']);
- grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma_{r,i}$')
-% ylim([0,2.0]);
- if(1)
- subplot(223)
- plot(ky,g_I(1,:),'-','DisplayName','Primar. instability'); hold on;
- plot(kx,g_II(:,1),'x-','DisplayName','Second. instability'); hold on;
- plot([max(ky)*2/3,max(ky)*2/3],[0,10],'--k', 'DisplayName','2/3 Orszag AA');
- grid on; xlabel('$k\rho_s$'); ylabel('$\gamma R/c_s$'); legend('show');
- ylim([0,max(g_I(1,:))]); xlim([0,max(ky)]);
- shearplot = 426; phiplot = 428;
- else
- shearplot = 223; phiplot = 224;
- end
- subplot(shearplot)
- clr = line_colors(min(ip,numel(line_colors(:,1))),:);
- lstyle = line_styles(min(ij,numel(line_styles)));
- plot(Ts2D,shear_maxr_maxz,'DisplayName','$\max_{r,z}(s)$'); hold on;
- plot(Ts2D,shear_maxr_avgz,'DisplayName','$\max_{r}\langle s \rangle_z$'); hold on;
- plot(Ts2D,shear_avgr_maxz,'DisplayName','$\max_{z}\langle s \rangle_r$'); hold on;
- plot(Ts2D,shear_avgr_avgz,'DisplayName','$\langle s \rangle_{r,z}$'); hold on;
- grid on; xlabel('$t c_s/R$'); ylabel('$shear$');
- subplot(phiplot)
- clr = line_colors(min(ip,numel(line_colors(:,1))),:);
- lstyle = line_styles(min(ij,numel(line_styles)));
- plot(Ts2D,phi_maxr_maxz,'DisplayName','$\max_{r,z}(\phi)$'); hold on;
- plot(Ts2D,phi_maxr_avgz,'DisplayName','$\max_{r}\langle\phi\rangle_z$'); hold on;
- plot(Ts2D,phi_avgr_maxz,'DisplayName','$\max_{z}\langle\phi\rangle_r$'); hold on;
- plot(Ts2D,phi_avgr_avgz,'DisplayName','$\langle\phi\rangle_{r,z}$'); hold on;
- grid on; xlabel('$t c_s/R$'); ylabel('$E.S. pot$');
-save_figure
-end
-
-if 1
-%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG = 5000; % Averaging time duration
-%Compute steady radial transport
-tend = Ts0D(end); tstart = tend - TAVG;
-[~,its0D] = min(abs(Ts0D-tstart));
-[~,ite0D] = min(abs(Ts0D-tend));
-SCALE = (2*pi/Nx/Ny)^2;
-gamma_infty_avg = mean(PGAMMA_RI(its0D:ite0D))*SCALE;
-gamma_infty_std = std (PGAMMA_RI(its0D:ite0D))*SCALE;
-% Compute steady shearing rate
-tend = Ts2D(end); tstart = tend - TAVG;
-[~,its2D] = min(abs(Ts2D-tstart));
-[~,ite2D] = min(abs(Ts2D-tend));
-shear_infty_avg = mean(shear_maxr_avgz(its2D:ite2D));
-shear_infty_std = std (shear_maxr_avgz(its2D:ite2D));
-Q_infty_avg = mean(Q_RI(its2D:ite2D))*SCALE;
-Q_infty_std = std(Q_RI(its2D:ite2D))*SCALE;
-% plots
-fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
- subplot(311)
-% yyaxis left
- plot(Ts0D,PGAMMA_RI*SCALE,'DisplayName','$\langle n_i d\phi/dz \rangle_z$'); hold on;
- plot(Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',...
- 'DisplayName',['$\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]);
- grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_r$')
- ylim([0,5*abs(gamma_infty_avg)]); xlim([0,Ts0D(end)]);
- title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
- ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
- ' $\mu_{hd}=$',num2str(MU)]);
-% yyaxis right
-% plot(Ts2D,Q_RI*SCALE,'.','DisplayName','$\langle T_i d\phi/dz \rangle_z$'); hold on;
-% ylim([0,5*Q_infty_avg]); xlim([0,Ts0D(end)]); ylabel('$Q_r$')
-% plot(Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Q_infty_avg, '--k',...
-% 'DisplayName',['$Q^{\infty} = $',num2str(Q_infty_avg),'$\pm$',num2str(Q_infty_std)]);
-% legend('show','Location','west')
- %
- subplot(312)
- clr = line_colors(1,:);
- lstyle = line_styles(1);
- plot(Ts2D,shear_maxr_maxz,'DisplayName','$\max_{r,z}(s_\phi)$'); hold on;
- plot(Ts2D,shear_maxr_avgz,'DisplayName','$\max_{r}\langle s_\phi\rangle_z$'); hold on;
- plot(Ts2D,shear_avgr_maxz,'DisplayName','$\max_{z}\langle s_\phi\rangle_r$'); hold on;
- plot(Ts2D,shear_avgr_avgz,'DisplayName','$\langle s_\phi\rangle_{r,z}$'); hold on;
- plot(Ts2D(its2D:ite2D),ones(ite2D-its2D+1,1)*shear_infty_avg, '-k',...
- 'DisplayName',['$s^{\infty} = $',num2str(shear_infty_avg),'$\pm$',num2str(shear_infty_std)]);
- ylim([0,shear_infty_avg*5.0]); xlim([0,Ts0D(end)]);
- grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show');
- subplot(313)
- [TY,TX] = meshgrid(r,Ts2D);
-% pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,:),2))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_r\phi\rangle_z$') %colorbar;
- pclr = pcolor(TX,TY,squeeze(mean(dr2phi(:,:,:),2))'); set(pclr, 'edgecolor','none'); legend('Shear ($\langle \partial_r^2\phi\rangle_z$)') %colorbar;
- caxis(1*shear_infty_avg*[-1 1]); xlabel('$t c_s/R$'), ylabel('$r/\rho_s$'); colormap(bluewhitered(256))
-save_figure
-end
-
-if 1
-%% Space time diagramms
-tstart = 0; tend = Ts2D(end);
-[~,itstart] = min(abs(Ts2D-tstart));
-[~,itend] = min(abs(Ts2D-tend));
-trange = itstart:itend;
-[TY,TX] = meshgrid(r,Ts2D(trange));
-fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
- subplot(211)
-% pclr = pcolor(TX,TY,squeeze(mean(dens_i(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
- pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
- shading interp
- colormap hot;
- caxis([0.0,0.05*max(max(mean(ni00(:,:,its2D:ite2D).*dzphi(:,:,its2D:ite2D),2)))]);
- caxis([0.0,0.01]); c = colorbar; c.Label.String ='\langle n_i\partial_z\phi\rangle_z';
- xticks([]); ylabel('$r/\rho_s$')
-% legend('Radial part. transport $\langle n_i\partial_z\phi\rangle_z$')
- title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
- ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
- ' $\mu_{hd}=$',num2str(MU)]);
-% subplot(312)
-% pclr = pcolor(TX,TY,squeeze(mean(temp_i(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
-% shading interp
-% % colormap(bluewhitered(256));
-% xticks([]); ylabel('$r/\rho_s$')
-% legend('Radial part. transport $\langle T_i\partial_z\phi\rangle_z$')
-% title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
-% ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
-% ' $\mu_{hd}=$',num2str(MU)]);
- subplot(212)
- pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
- fieldmax = max(max(mean(abs(drphi(:,:,its2D:ite2D)),2)));
- caxis([-fieldmax,fieldmax]); c = colorbar; c.Label.String ='\langle \partial_r\phi\rangle_z';
- xlabel('$t c_s/R$'), ylabel('$r/\rho_s$')
-% legend('Zonal flow $\langle \partial_r\phi\rangle_z$')
-save_figure
-end
-
-if 0
-%% Averaged shear and Reynold stress profiles
-trange = its2D:ite2D;
-% trange = 100:200;
-figure;
-plt = @(x) squeeze(mean(mean(real(x(:,:,trange)),2),3))/max(abs(squeeze(mean(mean(real(x(:,:,trange)),2),3))));
-plot(r,plt(-dr2phi),'-k','DisplayName','Zonal shear'); hold on;
-plot(r,plt(-drphi.*dzphi),'-','DisplayName','$\Pi_\phi$'); hold on;
-% plot(r,plt(-drphi.*dzTe),'-','DisplayName','$\Pi_{Te}$'); hold on;
-plot(r,plt(-drphi.*dzTi),'-','DisplayName','$\Pi_{Ti}$'); hold on;
-plot(r,plt(-drphi.*dzphi-drphi.*dzTi),'-','DisplayName','$\Pi_\phi+\Pi_{Ti}$'); hold on;
-% plot(r,plt(-drphi.*dzphi-drphi.*dzTi-drphi.*dzTe),'-','DisplayName','$\Pi_\phi+\Pi_{Te}+\Pi_{Ti}$'); hold on;
-xlim([-L/2,L/2]); xlabel('$r/\rho_s$'); grid on; legend('show')
-end
-
-if 1
-%% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
-tstart = 0.8*Ts2D(end); tend = Ts2D(end);
-[~,itstart] = min(abs(Ts2D-tstart));
-[~,itend] = min(abs(Ts2D-tend));
-trange = itstart:itend;
-%full kperp points
-phi_k_2 = reshape(mean((abs(PHI(:,:,trange))).^2,3),[numel(kx),1]);
-kperp = reshape(sqrt(kx.^2+ky.^2),[numel(kx),1]);
-% interpolated kperps
-nk_noAA = floor(2/3*numel(kx));
-kp_ip = kx;
-[thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip);
-[xn,yn] = pol2cart(thg,rg);
-[ky_s, sortIdx] = sort(ky);
-[xc,yc] = meshgrid(ky_s,kx);
-Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,trange))).^2,3)),xn,yn);
-phi_kp = mean(Z_rth,2);
-Z_rth = interp2(xc,yc,squeeze(mean((abs(Ni00(:,sortIdx,trange))).^2,3)),xn,yn);
-ni00_kp = mean(Z_rth,2);
-Z_rth = interp2(xc,yc,squeeze(mean((abs(Ne00(:,sortIdx,trange))).^2,3)),xn,yn);
-ne00_kp = mean(Z_rth,2);
-%for theorical trends
-a1 = phi_kp(2)*kp_ip(2).^(13/3);
-a2 = phi_kp(2)*kp_ip(2).^(3)./(1+kp_ip(2).^2).^(-2);
-fig = figure; FIGNAME = ['cascade','_',PARAMS];set(gcf, 'Position', [100, 100, 500, 500])
-% scatter(kperp,phi_k_2,'.k','MarkerEdgeAlpha',0.4,'DisplayName','$|\phi_k|^2$'); hold on; grid on;
-plot(kp_ip,phi_kp,'^','DisplayName','$\langle|\phi_k|^2\rangle_{k_\perp}$'); hold on;
-plot(kp_ip,ni00_kp, '^','DisplayName','$\langle|N_{i,k}^{00}|^2\rangle_{k_\perp}$'); hold on;
-plot(kp_ip,ne00_kp, '^','DisplayName','$\langle|N_{e,k}^{00}|^2\rangle_{k_\perp}$'); hold on;
-plot(kp_ip,a1*kp_ip.^(-13/3),'-','DisplayName','$k^{-13/3}$');
-plot(kp_ip,a2/100*kp_ip.^(-3)./(1+kp_ip.^2).^2,'-','DisplayName','$k^{-3}/(1+k^2)^2$');
-set(gca, 'XScale', 'log');set(gca, 'YScale', 'log'); grid on
-xlabel('$k_\perp \rho_s$'); legend('show','Location','southwest')
-title({['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
-', $L=',num2str(L),'$, $N=',num2str(Nx),'$'];[' $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
-' $\mu_{hd}=$',num2str(MU)]});
-xlim([0.1,10]);
-save_figure
-clear Z_rth phi_kp ni_kp Ti_kp
-end
-
-%%
-t0 =00;
-[~, it02D] = min(abs(Ts2D-t0));
-[~, it05D] = min(abs(Ts5D-t0));
-skip_ = 4;
-DELAY = 1e-3*skip_;
-FRAMES_2D = it02D:skip_:numel(Ts2D);
-FRAMES_5D = it05D:skip_:numel(Ts5D);
-%% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-if 0
-%% part density electron
-GIFNAME = ['ne',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-FIELD = real(dens_e); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
-FIELDNAME = '$n_e$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-create_gif
-% create_mov
-end
-if 0
-%% part temperature electron
-GIFNAME = ['Te',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-FIELD = real(temp_e); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
-FIELDNAME = '$T_e$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-create_gif
-% create_mov
-end
-if 0
-%% part density ion
-GIFNAME = ['ni',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-FIELD = real(dens_i); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
-FIELDNAME = '$n_i$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-create_gif
-% create_mov
-end
-if 0
-%% part temperature ion
-GIFNAME = ['Ti',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-FIELD = real(temp_i); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
-FIELDNAME = '$T_i$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-create_gif
-% create_mov
-end
-if 0
-%% GC Density ion
-GIFNAME = ['ni00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
-FIELDNAME = '$n_i^{00}$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-create_gif
-% create_mov
-end
-if 0
-%% GC Density electrons
-GIFNAME = ['ne00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-FIELD = real(ne00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
-FIELDNAME = '$n_e^{00}$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-% create_gif
-create_mov
-end
-if 0
-%% Phi real space
-GIFNAME = ['phi',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 0;
-FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
-FIELDNAME = '$\phi$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-create_gif
-% create_mov
-end
-if 0
-%% shear
-GIFNAME = ['shear_r',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 1;
-FIELD = -dr2phi; X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
-FIELDNAME = '$s$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-create_gif
-end
-if 0
-%% phi averaged on z
-GIFNAME = ['phi_z0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING=1;
-FIELD =(squeeze(mean(real(phi),2))); linestyle = '-.k'; FRAMES = FRAMES_2D;
-X = (r); T = Ts2D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r);
-FIELDNAME = '$\langle\phi\rangle_{z}$'; XNAME = '$r/\rho_s$';
-create_gif_1D_phi
-end
-if 0
-%% flow averaged on z
-GIFNAME = ['zf_z0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING=1;
-FIELD =(squeeze(mean(real(drphi),2))); linestyle = '-.k'; FRAMES = FRAMES_2D;
-X = (r); T = Ts2D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r);
-FIELDNAME = '$\langle\phi\rangle_{z}$'; XNAME = '$r/\rho_s$';
-create_gif_1D_phi
-end
-if 0
-%% shear averaged on z
-GIFNAME = ['shear_z0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING=1;
-FIELD =(squeeze(mean(real(dr2phi),2))); linestyle = '-.k'; FRAMES = FRAMES_2D;
-X = (r); T = Ts2D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r);
-FIELDNAME = '$\langle\phi\rangle_{z}$'; XNAME = '$r/\rho_s$';
-create_gif_1D_phi
-end
-if 0
-%% phi kperp spectrum
-GIFNAME = ['phi2_kp',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING = 0;
-FIELD =log10(phi_kp_t); linestyle = '-'; FRAMES = FRAMES_2D;
-X = kp_ip; T = Ts2D; YMIN = -20; YMAX = 10; XMIN = min(kx); XMAX = max(kx);
-FIELDNAME = '$\log_{10}|\tilde\phi_k|^2$'; XNAME = '$k_r\rho_s$';
-create_gif_1D
-end
-if 0
-%% Density ion frequency
-GIFNAME = ['Ni00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; FRAMES = FRAMES_2D;
-FIELD =ifftshift((abs(Ni00)),2); X = fftshift(kx,2); Y = fftshift(ky,2); T = Ts2D;
-FIELDNAME = '$N_i^{00}$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$';
-create_gif
-end
-if 0
-%% Density electron frequency
-GIFNAME = ['Ne00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; FRAMES = FRAMES_2D;
-FIELD =ifftshift((abs(Ne00)),2); X = fftshift(kx,2); Y = fftshift(ky,2); T = Ts2D;
-FIELDNAME = '$N_e^{00}$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$';
-create_gif
-end
-if 0
-%% kx vs P Si
-GIFNAME = ['Sip0_kx',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-plt = @(x) squeeze(max((abs(x)),[],4));
-FIELD =plt(Sipj(:,1,:,:,:)); X = kx'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
-FIELDNAME = '$N_i^{p0}$'; XNAME = '$k_{max}\rho_s$'; YNAME = '$P$, $\sum_z$';
-create_gif_imagesc
-end
-if 0
-%% maxky, kx vs p, for all Nipj over time
-GIFNAME = ['Nipj_kx',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-plt = @(x) squeeze(sum((abs(x)),4));
-FIELD = plt(Nipj); X = kx'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
-FIELDNAME = 'N_i'; XNAME = '$k_r\rho_s$'; YNAME = '$P$';
-create_gif_5D
-end
-if 0
-%% maxkx, ky vs p, for all Nipj over time
-GIFNAME = ['Nipj_ky',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-plt = @(x) fftshift(squeeze(sum((abs(x)),3)),3);
-FIELD = plt(Nipj); X = sort(ky'); Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
-FIELDNAME = 'N_i'; XNAME = '$k_z\rho_s$'; YNAME = '$P$, $\sum_r$';
-create_gif_5D
-end
-%%
-ZF_fourier_analysis
-end
\ No newline at end of file
diff --git a/wk/analysis_stat_2D.m b/wk/analysis_stat_2D.m
deleted file mode 100644
index cae39cf8a58e7212bdf2e5c0b615dd4407684a52..0000000000000000000000000000000000000000
--- a/wk/analysis_stat_2D.m
+++ /dev/null
@@ -1,53 +0,0 @@
-if 0
-%% ni ne phase space for different time windows
-fig = figure;
- t0 = 3500; t1 = 4500; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
- plt = @(x) abs(squeeze(reshape(x(2,:,it0:it1),[],1)));
- plot(plt(Ni00),plt(Ne00),'.');
- title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
- xlabel('$n_i$'); ylabel('$n_e$');
-end
-
-if 0
-%% histograms
-t0 = 2700; t1 = 2800; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
-plt = @(x) squeeze(reshape(x(:,:,it0:it1),[],1));
-fig = figure('Position', [100, 100, 800, 500]);
-subplot(221)
- hist(plt(ne00),50); hold on
- title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
- xlabel('$n_e$');
-subplot(222)
- hist(plt(ni00),50); hold on
- title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
- xlabel('$n_i$');
-subplot(223)
- hist(plt(phi),50); hold on
- title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
- xlabel('$\phi$');
-FIGNAME = ['histograms_ne_ni_phi'];
-save_figure
-end
-
-if 0
-%% ni ne phi 3D phase space
-fig = figure;
- t0 = 2700; [~,it0] = min(abs(t0-Ts2D));
- plt = @(x) squeeze(reshape(x(:,:,it0),[],1));
- plot3(plt(ni00),plt(ne00),plt(phi),'.','MarkerSize',1.9);
- title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
- xlabel('$n_i$'); ylabel('$n_e$'); zlabel('$\phi$'); grid on;
-end
-
-t0 = 2400;
-skip_ = 2;
-DELAY = 0.01*skip_;
-FRAMES = floor(t0/(Ts2D(2)-Ts2D(1)))+1:skip_:numel(Ts2D);
-if 0
-%% ni ne phi 3D phase space GIF
-GIFNAME = ['ni_ne_phi_phase_space',sprintf('_%.2d',JOBNUM)];
-X = ne00; Y = ni00; Z = phi; T = Ts2D;
-MARKERSIZE = 0.01;
-XNAME = '$n_i$'; YNAME = '$n_e$'; ZNAME = '$\phi$'; VIEW = [1,-1,1];
-create_gif_3D
-end
\ No newline at end of file
diff --git a/wk/compile_cosolver_mat.m b/wk/compile_cosolver_mat.m
deleted file mode 100644
index 7b2310dd0556e38a05b6b6a1f155383d9735d95c..0000000000000000000000000000000000000000
--- a/wk/compile_cosolver_mat.m
+++ /dev/null
@@ -1,137 +0,0 @@
-addpath(genpath('../matlab')) % ... add
-%% Grid configuration
-N = 200; % Frequency gridpoints (Nkx = N/2)
-L = 120; % Size of the squared frequency domain
-dk = 2.*pi/L;
-kmax = N/2*dk;
-kx = dk*(0:N/2);
-ky = dk*(0:N/2);
-[ky, kx]= meshgrid(ky,kx);
-KPERP = sqrt(kx.^2 + ky.^2);
-kperp = reshape(KPERP,[1,numel(kx)^2]);
-kperp = uniquetol(kperp,1e-14);
-Nperp = numel(kperp);
-%% Model
- % ! 0 = nothing
- % ! 1 = coulomb
- % ! 2 = pitch-angle (with constant coll.freq.)
- % ! 3 = sugama
- % ! 4 = pitch-angle with veloty dependent coll. freq.
- % ! 5 = improved sugama
- % ! 6 = Hirschman-Sigmar Clarke
- % ! 7 = Abel (only for like species)
- % ! 8 = conserving momentun pitch-angle operator (with veloty dependent coll. freq.)
- % ! 9 = GK coulomb polarization term
-CO = 3;
-GK = 1;
-P = 10;
-J = 5;
-M_FLR= 0; %to increase the NFLR
-if 0
- %% plot the kperp distribution
- figure
- plot(kperp)
-end
-% %% Check if the differences btw kperp is larger than naming precision
-%%
-%% We compute only on a kperp grid with dk space from 0 to kperpmax
-kperpmax = sqrt(2) * kmax;
-kperp = unique([0:dk:kperpmax,kperpmax]);
-%% Naming
-if CO == 1; CONAME = 'FC'; end;
-if CO == 2; CONAME = 'PA'; end;
-if CO == 3; CONAME = 'SG'; end;
-% matfilename = ['../iCa/',CONAME,'_P_',num2str(P),'_J_',num2str(J),...
-% '_N_',num2str(N),'_dk_',num2str(dk),'_MFLR_',num2str(M_FLR),'.h5'];
-matfilename = '../iCa/gk_sugama_P_10_J_5.h5';
-n_ = 0;
-for k_ = kperp
- disp(['-Writing matrices for kperp = ',num2str(k_)])
- %% Script to run COSOlver in order to create needed collision matrices
- COSOlver_path = '../../Documents/MoliSolver/COSOlver/';
- COSOLVER.pmaxe = P;
- COSOLVER.jmaxe = J;
- COSOLVER.pmaxi = P;
- COSOLVER.jmaxi = J;
- COSOLVER.kperp = k_;
-
- COSOLVER.neFLR = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2)))+M_FLR; % rule of thumb for sum truncation
- COSOLVER.niFLR = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2)))+M_FLR;
- COSOLVER.idxT4max = 40;
-
- COSOLVER.neFLRs = 0; % ... only for GK abel
- COSOLVER.npeFLR = 0; % ... only for GK abel
- COSOLVER.niFLRs = 0; % ... only for GK abel
- COSOLVER.npiFLR = 0; % ... only for GK abel
-
- COSOLVER.gk = GK;
- COSOLVER.co = CO;
-
- k_string = sprintf('%0.4f',k_);
- n_string = sprintf('%0.5d',n_);
- self_mat_file_name = ['../iCa/self_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
- '_ESELF_',num2str(COSOLVER.co),'_ISELF_',num2str(COSOLVER.co),...
- '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
- '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
- '_JE_12',...
- '_NFLR_',num2str(COSOLVER.neFLR),...
- '_kperp_',k_string,'.h5'];
- ei_mat_file_name = ['../iCa/ei_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
- '_ETEST_',num2str(COSOLVER.co),'_EBACK_',num2str(COSOLVER.co),...
- '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
- '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
- '_JE_12_tau_1.0000_mu_0.0233',...
- '_NFLRe_',num2str(COSOLVER.neFLR),'_NFLRi_',num2str(COSOLVER.neFLR)...
- '_kperp_',k_string,'.h5'];
- ie_mat_file_name = ['../iCa/ie_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
- '_ITEST_',num2str(COSOLVER.co),'_IBACK_',num2str(COSOLVER.co),...
- '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
- '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
- '_JE_12_tau_1.0000_mu_0.0233',...
- '_NFLRe_',num2str(COSOLVER.neFLR),'_NFLRi_',num2str(COSOLVER.neFLR)...
- '_kperp_',k_string,'.h5'];
-
- %% Load self matrix
- C_self = h5read(self_mat_file_name,'/Caapj/Ceepj');
- sz_ = size(C_self);
- % Write it in the compiled file
- h5create(matfilename,['/',n_string,'/Caapj/Ceepj'],sz_)
- h5write (matfilename,['/',n_string,'/Caapj/Ceepj'],C_self)
- h5create(matfilename,['/',n_string,'/Caapj/Ciipj'],sz_)
- h5write (matfilename,['/',n_string,'/Caapj/Ciipj'],C_self)
- %% Load ei matrices
- % Field
- C_eiF = h5read(ei_mat_file_name,'/Ceipj/CeipjF');
- sz_ = size(C_eiF);
- h5create(matfilename,['/',n_string,'/Ceipj/CeipjF'],sz_)
- h5write (matfilename,['/',n_string,'/Ceipj/CeipjF'],C_eiF)
- % Test
- C_eiT = h5read(ei_mat_file_name,'/Ceipj/CeipjT');
- sz_ = size(C_eiT);
- h5create(matfilename,['/',n_string,'/Ceipj/CeipjT'],sz_)
- h5write (matfilename,['/',n_string,'/Ceipj/CeipjT'],C_eiT)
-
- %% Load ie matrices
- % Field
- C_ieF = h5read(ie_mat_file_name,'/Ciepj/CiepjF');
- sz_ = size(C_ieF);
- h5create(matfilename,['/',n_string,'/Ciepj/CiepjF'],sz_)
- h5write (matfilename,['/',n_string,'/Ciepj/CiepjF'],C_ieF)
- % Test
- C_ieT = h5read(ie_mat_file_name,'/Ciepj/CiepjT');
- sz_ = size(C_eiT);
- h5create(matfilename,['/',n_string,'/Ciepj/CiepjT'],sz_)
- h5write (matfilename,['/',n_string,'/Ciepj/CiepjT'],C_ieT)
-
- %% Copy fort.90 input file and put grid params
- if(k_ == 0)
- h5create(matfilename,'/coordkperp',numel(kperp));
- h5write (matfilename,'/coordkperp',kperp);
- h5create(matfilename,'/dims_e',2);
- h5write (matfilename,'/dims_e',[P,J]);
- h5create(matfilename,'/dims_i',2);
- h5write (matfilename,'/dims_i',[P,J]);
- end
- n_ = n_ + 1;
-end
-disp(['File saved @ :',matfilename])
\ No newline at end of file
diff --git a/wk/compute_collision_mat.m b/wk/compute_collision_mat.m
deleted file mode 100644
index ac5214c15838a5ea489ab1a71fb08f016696d6b5..0000000000000000000000000000000000000000
--- a/wk/compute_collision_mat.m
+++ /dev/null
@@ -1,99 +0,0 @@
-addpath(genpath('../matlab')) % ... add
-%% Grid configuration
-N = 200; % Frequency gridpoints (Nkx = N/2)
-L = 120; % Size of the squared frequency domain
-dk = 2.*pi/L;
-kmax = N/2*dk;
-kx = dk*(0:N/2);
-ky = dk*(0:N/2);
-[ky, kx]= meshgrid(ky,kx);
-KPERP = sqrt(kx.^2 + ky.^2);
-kperp = reshape(KPERP,[1,numel(kx)^2]);
-kperp = uniquetol(kperp,1e-14);
-Nperp = numel(kperp);
-%% Model
- % ! 0 = nothing
- % ! 1 = coulomb
- % ! 2 = pitch-angle (with constant coll.freq.)
- % ! 3 = sugama
- % ! 4 = pitch-angle with veloty dependent coll. freq.
- % ! 5 = improved sugama
- % ! 6 = Hirschman-Sigmar Clarke
- % ! 7 = Abel (only for like species)
- % ! 8 = conserving momentun pitch-angle operator (with veloty dependent coll. freq.)
- % ! 9 = GK coulomb polarization term
-CO = 3;
-GK = 1;
-
-%%
-%% We compute only on a kperp grid with dk space from 0 to kperpmax
-kperpmax = sqrt(2) * kmax;
-kperp = unique([0:dk:kperpmax,kperpmax]);
-%%
-%% If DK collision, compute only kperp = 0
-if (GK == 0); kperp = 0; end;
-%% Kperp scan
-n_ = 1;
-for k_ = kperp
- disp(['----------Computing matrix ',num2str(n_),'/',num2str(numel(kperp)),'----------'])
- %% Script to run COSOlver in order to create needed collision matrices
- COSOlver_path = '../../Documents/MoliSolver/COSOlver/';
- COSOLVER.pmaxe = 10;
- COSOLVER.jmaxe = 05;
- COSOLVER.pmaxi = 10;
- COSOLVER.jmaxi = 05;
- COSOLVER.kperp = k_;
-
- COSOLVER.neFLR = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2))); % rule of thumb for sum truncation
- COSOLVER.niFLR = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2)));
- COSOLVER.idxT4max = 40;
-
- COSOLVER.neFLRs = 0; % ... only for GK abel
- COSOLVER.npeFLR = 0; % ... only for GK abel
- COSOLVER.niFLRs = 0; % ... only for GK abel
- COSOLVER.npiFLR = 0; % ... only for GK abel
-
- COSOLVER.gk = GK;
- COSOLVER.co = CO;
-
- write_fort90_COSOlver
-
- k_string = sprintf('%0.4f',k_);
- self_mat_file_name = ['../iCa/self_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
- '_ESELF_',num2str(COSOLVER.co),'_ISELF_',num2str(COSOLVER.co),...
- '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
- '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
- '_JE_12',...
- '_NFLR_',num2str(COSOLVER.neFLR),...
- '_kperp_',k_string,'.h5'];
- ei_mat_file_name = ['../iCa/ei_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
- '_ETEST_',num2str(COSOLVER.co),'_EBACK_',num2str(COSOLVER.co),...
- '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
- '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
- '_JE_12_tau_1.0000_mu_0.0233',...
- '_NFLRe_',num2str(COSOLVER.neFLR),'_NFLRi_',num2str(COSOLVER.neFLR)...
- '_kperp_',k_string,'.h5'];
- ie_mat_file_name = ['../iCa/ie_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
- '_ITEST_',num2str(COSOLVER.co),'_IBACK_',num2str(COSOLVER.co),...
- '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
- '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
- '_JE_12_tau_1.0000_mu_0.0233',...
- '_NFLRe_',num2str(COSOLVER.neFLR),'_NFLRi_',num2str(COSOLVER.neFLR)...
- '_kperp_',k_string,'.h5'];
-% if (exist(self_mat_file_name,'file')>0 && exist(ei_mat_file_name,'file')>0 && exist(ie_mat_file_name,'file') > 0)
- if (exist(ei_mat_file_name,'file')>0)%&& exist(ie_mat_file_name,'file') > 0)
- disp(['Matrix available for kperp = ',k_string]);
- else
- cd ../../Documents/MoliSolver/COSOlver/
- disp(['Matrix not found for kperp = ',k_string]);
- disp([num2str(n_),'/',Nperp])
- disp('computing...');
- CMD = 'mpirun -np 6 bin/CO 2 2 2 > out.txt';
- disp(CMD);
- system(CMD);
- system(CMD);
- disp('..done');
- cd ../../../HeLaZ/wk
- end
- n_ = n_ + 1;
-end
\ No newline at end of file
diff --git a/wk/daint_run.m b/wk/daint_run.m
deleted file mode 100644
index 95806c96bdeb2983e39b1a04085193192006852e..0000000000000000000000000000000000000000
--- a/wk/daint_run.m
+++ /dev/null
@@ -1,89 +0,0 @@
-%clear all;
-addpath(genpath('../matlab')) % ... add
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Set Up parameters
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% CLUSTER PARAMETERS
-CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss
-NP_P = 2; % MPI processes along p
-NP_KX = 18; % MPI processes along kx
-CLUSTER.PART = 'normal'; % debug or normal
-if(strcmp(CLUSTER.PART,'debug')); CLUSTER.TIME = '00:30:00'; end;
-CLUSTER.MEM = '12GB'; % Memory
-CLUSTER.JNAME = 'gamma_inf';% Job name
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
-ETAB = 0.6; % Magnetic gradient
-NU_HYP = 1.0; % Hyperdiffusivity coefficient
-%% GRID PARAMETERS
-N = 200; % Frequency gridpoints (Nkx = N/2)
-L = 120; % Size of the squared frequency domain
-P = 12; % Electron and Ion highest Hermite polynomial degree
-J = 06; % Electron and Ion highest Laguerre polynomial degree
-MU_P = 0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J = 0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
-%% TIME PARAMETERS
-TMAX = 1000; % Maximal time unit
-DT = 1e-2; % Time step
-SPS0D = 1; % Sampling per time unit for profiler
-SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 1/40; % Sampling per time unit for 5D arrays
-SPSCP = 0; % Sampling per time unit for checkpoints
-RESTART = 1; % To restart from last checkpoint
-JOB2LOAD= 2;
-%% OPTIONS
-SIMID = ['HeLaZ_v2.5_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation
-% SIMID = 'test_marconi_sugama'; % Name of the simulation
-SIMID = sprintf(SIMID,NU);
-PREFIX =[];
-% PREFIX = sprintf('%d_%d_',NP_P, NP_KX);
-% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
-CO = 1;
-CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
-NL_CLOS = 1; % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS)
-KERN = 0; % Kernel model (0 : GK)
-INIT_PHI= 1; % Start simulation with a noisy phi and moments
-%% OUTPUTS
-W_DOUBLE = 1;
-W_GAMMA = 1;
-W_PHI = 1;
-W_NA00 = 1;
-W_NAPJ = 1;
-W_SAPJ = 0;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% fixed parameters (for current study)
-KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
-KXEQ0 = 0; % put kx = 0
-KPAR = 0.0; % Parellel wave vector component
-LAMBDAD = 0.0;
-NON_LIN = 1 *(1-KXEQ0); % activate non-linearity (is cancelled if KXEQ0 = 1)
-PMAXE = P; % Highest electron Hermite polynomial degree
-JMAXE = J; % Highest '' Laguerre ''
-PMAXI = P; % Highest ion Hermite polynomial degree
-JMAXI = J; % Highest '' Laguerre ''
-kmax = N*pi/L;% Highest fourier mode
-HD_CO = 0.5; % Hyper diffusivity cutoff ratio
-MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
-NOISE0 = 1.0e-5;
-ETAT = 0.0; % Temperature gradient
-ETAN = 1.0; % Density gradient
-TAU = 1.0; % e/i temperature ratio
-% Compute processes distribution
-Ntot = NP_P * NP_KX;
-Nnodes = ceil(Ntot/36);
-Nppn = Ntot/Nnodes;
-CLUSTER.NODES = num2str(Nnodes); % MPI process along p
-CLUSTER.NTPN = num2str(Nppn); % MPI process along kx
-CLUSTER.CPUPT = '1'; % CPU per task
-CLUSTER.NTPC = '1'; % N tasks per core (openmp threads)
-%% Run file management scripts
-setup
-write_sbash_daint
-system('rm fort.90 setup_and_run.sh batch_script.sh');
-disp('done');
-if(mod(NP_P*NP_KX,36)~= 0)
- disp('WARNING : unused cores (ntot cores must be a 36 multiple)');
-end
\ No newline at end of file
diff --git a/wk/izar_run.m b/wk/izar_run.m
deleted file mode 100644
index 035adff386eabe7702c1f1d0c5c6ba7f28c4ebdf..0000000000000000000000000000000000000000
--- a/wk/izar_run.m
+++ /dev/null
@@ -1,86 +0,0 @@
-%clear all;
-addpath(genpath('../matlab')) % ... add
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Set Up parameters
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% CLUSTER PARAMETERS
-CLUSTER.TIME = '00:15:00'; % allocation time hh:mm:ss
-CLUSTER.PART = 'gpu'; % debug/gpu
-CLUSTER.MEM = '4G'; % Memory
-CLUSTER.JNAME = 'gamma_inf';% Job name
-NP_P = 1; % MPI processes along p (nodes)
-NP_KX = 20; % MPI processes along kx (cpu)
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
-ETAB = 0.6; % Magnetic gradient
-NU_HYP = 1.0; % Hyperdiffusivity coefficient
-%% GRID PARAMETERS
-N = 50; % Frequency gridpoints (Nkx = N/2)
-L = 10; % Size of the squared frequency domain
-P = 04; % Electron and Ion highest Hermite polynomial degree
-J = 04; % Electron and Ion highest Laguerre polynomial degree
-MU_P = 0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J = 0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
-%% TIME PARAMETERS
-TMAX = 200; % Maximal time unit
-DT = 1e-2; % Time step
-SPS0D = 1; % Sampling per time unit for profiler
-SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 1/50; % Sampling per time unit for 5D arrays
-SPSCP = 0; % Sampling per time unit for checkpoints
-RESTART = 1; % To restart from last checkpoint
-JOB2LOAD= 0;
-%% OPTIONS
-% SIMID = ['HeLaZ_v2.5_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation
-% SIMID = sprintf(SIMID,NU);
-SIMID = 'izar_setup'; % Name of the simulation
-PREFIX =[];
-% PREFIX = sprintf('%d_%d_',NP_P, NP_KX);
-% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
-CO = 1;
-CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
-NL_CLOS = 1; % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS)
-KERN = 0; % Kernel model (0 : GK)
-INIT_PHI= 1; % Start simulation with a noisy phi and moments
-%% OUTPUTS
-W_DOUBLE = 1;
-W_GAMMA = 1;
-W_PHI = 1;
-W_NA00 = 1;
-W_NAPJ = 1;
-W_SAPJ = 0;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% fixed parameters (for current study)
-KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
-KXEQ0 = 0; % put kx = 0
-KPAR = 0.0; % Parellel wave vector component
-LAMBDAD = 0.0;
-NON_LIN = 1 *(1-KXEQ0); % activate non-linearity (is cancelled if KXEQ0 = 1)
-PMAXE = P; % Highest electron Hermite polynomial degree
-JMAXE = J; % Highest '' Laguerre ''
-PMAXI = P; % Highest ion Hermite polynomial degree
-JMAXI = J; % Highest '' Laguerre ''
-kmax = N*pi/L;% Highest fourier mode
-HD_CO = 0.5; % Hyper diffusivity cutoff ratio
-MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
-NOISE0 = 1.0e-5;
-ETAT = 0.0; % Temperature gradient
-ETAN = 1.0; % Density gradient
-TAU = 1.0; % e/i temperature ratio
-% Compute processes distribution
-Ntot = NP_P * NP_KX;
-Nnodes = ceil(Ntot/48);
-Nppn = Ntot/Nnodes;
-CLUSTER.NODES = num2str(Nnodes); % MPI process along p
-CLUSTER.NTPN = num2str(Nppn); % MPI process along kx
-CLUSTER.CPUPT = '1'; % CPU per task
-%% Run file management scripts
-setup
-write_sbash_izar
-system('rm fort.90 setup_and_run.sh batch_script.sh');
-disp('done');
-if(mod(NP_P*NP_KX,20)~= 0)
- disp('WARNING : unused cores (ntot cores must be a 20 multiple)');
-end
\ No newline at end of file
diff --git a/wk/ppb110_run.m b/wk/ppb110_run.m
deleted file mode 100644
index 2eddbcd16dbc25e17408aea3333b339b61913a95..0000000000000000000000000000000000000000
--- a/wk/ppb110_run.m
+++ /dev/null
@@ -1,106 +0,0 @@
-%clear all;
-addpath(genpath('../matlab')) % ... add
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Set Up parameters
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% CLUSTER PARAMETERS
-CLUSTER.TIME = '01:00:00'; % allocation time hh:mm:ss
-CLUSTER.NODES = '1'; % MPI process
-CLUSTER.CPUPT = '1'; % CPU per task
-CLUSTER.NTPN = '4'; % N tasks per node
-CLUSTER.MEM = '2000'; % Memory
-CLUSTER.JNAME = 'test_HeLaZ'; % Job name
-USERNAME = 'ahoffman'; % username at ppb110 for folder naming
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% PHYSICAL PARAMETERS
-NU = 1.0; % Collision frequency
-ETAB = 0.5; % Magnetic gradient
-NU_HYP = 0.1; % Hyperdiffusivity coefficient
-%% GRID PARAMETERS
-N = 150; % Frequency gridpoints (Nkx = N/2)
-L = 70; % Size of the squared frequency domain
-P = 2; % Electron and Ion highest Hermite polynomial degree
-J = 1; % Electron and Ion highest Laguerre polynomial degree
-MU_P = 0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J = 0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
-%% TIME PARAMETERS
-TMAX = 150; % Maximal time unit
-DT = 2e-2; % Time step
-SPS0D = 1; % Sampling per time unit for profiler
-SPS2D = 1/10; % Sampling per time unit for 2D arrays
-SPS5D = 1/10; % Sampling per time unit for 5D arrays
-SPSCP = 0; % Sampling per time unit for checkpoints
-RESTART = 0; % To restart from last checkpoint
-JOB2LOAD= 0;
-%% OPTIONS
-SIMID = 'test_low_sampling'; % Name of the simulation
-SIMID = sprintf(SIMID,NU);
-CO = -3; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
-CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
-KERN = 0; % Kernel model (0 : GK)
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% fixed parameters (for current study)
-KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. (not implemented)
-KXEQ0 = 0; % put kx = 0
-KPAR = 0.0; % Parellel wave vector component
-LAMBDAD = 0.0;
-NON_LIN = 1 *(1-KXEQ0); % activate non-linearity (is cancelled if KXEQ0 = 1)
-PMAXE = P; % Highest electron Hermite polynomial degree
-JMAXE = J; % Highest '' Laguerre ''
-PMAXI = P; % Highest ion Hermite polynomial degree
-JMAXI = J; % Highest '' Laguerre ''
-kmax = N*pi/L;% Highest fourier mode
-HD_CO = 0.5; % Hyper diffusivity cutoff ratio
-MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
-NOISE0 = 1.0e-5;
-ETAT = 0.0; % Temperature gradient
-ETAN = 1.0; % Density gradient
-TAU = 1.0; % e/i temperature ratio
-%% Run file management scripts
-setup % Write the input script "fort.90" with desired parameters
-
-%% Write the sh script to launch the job on slurm for PPB110
-INPUT = 'setup_and_run.sh';
-fid = fopen(INPUT,'wt');
-SCRATCH_SIMDIR = ['/scratch/',USERNAME,'/HeLaZ']; % Path to the simulation directory in the scratch
-% Writing of the script
-fprintf(fid,[...
-'#!/bin/bash\n',...
-'mkdir -p ',SCRATCH_SIMDIR,'/wk\n',...
-...
-'cd ',SCRATCH_SIMDIR,'/wk\n',...
-...
-'mkdir -p ', BASIC.RESDIR,'\n',...
-'cd ',BASIC.RESDIR,'\n',...
-'cp $HOME/HeLaZ/wk/fort.90 .\n',...
-'cp $HOME/HeLaZ/wk/batch_script.sh .\n',...
-...
-'sbatch batch_script.sh\n',...
-'echo $',SCRATCH_SIMDIR,'/results/',BASIC.SIMID,'/',BASIC.PARAMS,'/out.txt']);
-
-fclose(fid);
-system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']);
-
-%% Write the sbatch script for PPB110
-INPUT = 'batch_script.sh';
-fid = fopen(INPUT,'wt');
-
-fprintf(fid,[...
-'#!/bin/bash\n',...
-'#SBATCH --job-name=',CLUSTER.JNAME,'\n',...
-'#SBATCH --time=', CLUSTER.TIME,'\n',...
-'#SBATCH --nodes=', CLUSTER.NODES,'\n',...
-'#SBATCH --cpus-per-task=', CLUSTER.CPUPT,'\n',...
-'#SBATCH --ntasks-per-node=', CLUSTER.NTPN,'\n',...
-'#SBATCH --mem-per-cpu=', CLUSTER.MEM,'\n',...
-'#SBATCH --error=err.txt\n',...
-'#SBATCH --output=out.txt\n',...
-'module purge\n',...
-'module load PrgEnv-intel/17.0\n',...
-'srun ./../../../bin/helaz']);
-
-fclose(fid);
-system(['cp batch_script.sh ',BASIC.RESDIR,'/.']);
-
-disp('done');
diff --git a/wk/scaling_results.m b/wk/scaling_results.m
deleted file mode 100644
index 1ad96fa61aedff5d11c5219dd668749cad167246..0000000000000000000000000000000000000000
--- a/wk/scaling_results.m
+++ /dev/null
@@ -1,104 +0,0 @@
-default_plots_options
-
-% NPS = [01 02 04 08 12 16 20 24];
-NPS = [02 04 10];
-TIMES = NPS;
-if 0
-%% Load times
-i_ = 1;
-nn = 2;
-for np = NPS
- SIM_NAME = sprintf('nn%02d_np%02d',nn,np)
- hostfile = ['/marconi_scratch/userexternal/ahoffman/HeLaZ/results/',SIM_NAME,'/',BASIC.PARAMS];
- localtarget= ['../results/',sprintf('Scaling/nn%02d_np%02d',nn,np),'/'];
- system(['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',localtarget]);
- filename = ['../results/',sprintf('Scaling/nn%02d_np%02d',nn,np),'/',BASIC.PARAMS,'/outputs_00.h5'];
- TIMES(i_) = h5readatt(filename,'/data/input','cpu_time');
- i_ = i_ + 1;
-end
-disp(PARAMS)
-disp(TIMES')
-end
-%% Strong scaling measurement (no diagnostics)
-
-% Handwritten results for 512x256, P,J=2,1, Tmax = 10, mu=0, dt = 5e-2
-Results_512_21.np = NPS;
-Results_512_21.tproc = [857 447 217 107 74 59 46 41]; %Nthreads
-Results_512_21.tnode = [854 442 199 96 75 64 57 55]; %Nnodes
-
-% Handwritten results for 512x256, P,J=3,2, Tmax = 10, mu=0, dt 5e-2
-Results_512_32.np = [ 1, 2, 4, 8, 12, 16, 20, 24];
-Results_512_32.tproc = [2546 1290 0630 0316 0220 0177 0143 0125];
-
-% Handwritten results for 512x256, P,J=6,4, Tmax = 5, mu=0, dt 5e-2
-Results_512_64.np = [ 1, 2, 4, 8, 12, 16, 20, 24];
-Results_512_64.tproc = [5801 2974 1467 0746 0507 0408 0320 0279];
-
-% Handwritten results for 1024x512, P,J=2,1, Tmax = 5 dt = 0.05, mu = 0
-Results_1024_21.np = [ 1, 2, 4, 8, 12, 16, 20, 24];
-Results_1024_21.tproc = [2051, 1296, 0583, 0296, 0205, 0163, 0132, 0115];
-Results_1024_21.tnode = [2052 1327 0511 0231 0157 0109 0099 0083]; %Nnodes
-
-% Handwritten results for 1024x512, P,J=3,2, Tmax = 2 dt = 0.05, mu = 0
-Results_1024_32.np = [ 1, 2, 4, 8, 12, 16, 20, 24];
-Results_1024_32.tproc = [2248, 1526, 0644, 0326, 0233, 0182, 0148, 0130];
-
-% Handwritten results for 1024x512, P,J=6,4, Tmax = 2 dt = 0.05, mu = 0
-Results_1024_64.np = [ 1, 2, 4, 8, 12, 16, 20, 24];
-Results_1024_64.tproc = [10330,6548, 2901, 1497, 1001, 0769, 0610, 0520];
-
-%
-fig = figure;
-
- subplot(121)
- plot(1:24,1:24,'-k','DisplayName','Ideal')
- hold on
- res = Results_512_21;
- plot(res.np,res.tproc(1)./(res.tproc),'v-','Color',line_colors(1,:),...
- 'DisplayName','$512\times256$, $P,J=2,1$');
-
- res = Results_512_32;
- plot(res.np,res.tproc(1)./(res.tproc),'v-','Color',line_colors(2,:),...
- 'DisplayName','$512\times256$, $P,J=3,2$');
-
- res = Results_512_64;
- plot(res.np,res.tproc(1)./(res.tproc),'v-','Color',line_colors(3,:),...
- 'DisplayName','$512\times256$, $P,J=6,4$');
-
- res = Results_1024_21;
- plot(res.np,res.tproc(1)./(res.tproc),'^--','Color',line_colors(1,:),...
- 'DisplayName','$1024\times512$, $P,J=2,1$');
-
- res = Results_1024_32;
- plot(res.np,res.tproc(1)./(res.tproc),'^--','Color',line_colors(2,:),...
- 'DisplayName','$1024\times512$, $P,J=3,2$');
-
- res = Results_1024_64;
- plot(res.np,res.tproc(1)./(res.tproc),'^--','Color',line_colors(3,:),...
- 'DisplayName','$1024\times512$, $P,J=6,4$');
-
- xlabel('$N_p$'); ylabel('speedup')
- xlim([1,24]);
- legend('show')
- title('$N_n=01$')
-
-subplot(122)
- plot(1:24,1:24,'-k','DisplayName','Ideal')
- hold on
- res = Results_512_21;
- plot(res.np,res.tnode(1)./(res.tnode),'o-','Color',line_colors(1,:),...
- 'DisplayName','$512\times256$, $P,J=2,1$');
-
- res = Results_1024_21;
- plot(res.np,res.tnode(1)./(res.tnode),'o--','Color',line_colors(1,:),...
- 'DisplayName','$1024\times512$, $P,J=2,1$');
-
- xlabel('$N_n$'); ylabel('speedup')
- xlim([1,24]);
- legend('show')
- title('$N_p=01$')
-
-
-FIGNAME = '/home/ahoffman/HeLaZ/results/strong_scaling_new.pdf';
-saveas(fig,FIGNAME);
-disp(['Figure saved @ : ',FIGNAME])
diff --git a/wk/spectral_analysis.m b/wk/spectral_analysis.m
deleted file mode 100644
index 8830f118d3566211c8096cf531fa8a33cdae2819..0000000000000000000000000000000000000000
--- a/wk/spectral_analysis.m
+++ /dev/null
@@ -1,67 +0,0 @@
-
-%%
-if 0
-%% Hermite energy spectra
-% tf = Ts2D(end-3);
-skip = 1;
-fig = figure; FIGNAME = ['hermite_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 1800, 600]);
-plt = @(x) squeeze(x);
-N10 = floor(Nji/10);
-Nrow = (1+N10);
-Ncol = ceil(Nji/Nrow);
-for ij = 1:Nji
- subplot(Nrow,Ncol,ij)
- for it5 = 1:skip:Ns5D
- alpha = it5*1.0/Ns5D;
- loglog(Pi(1:2:end),plt(epsilon_i_pj(1:2:end,ij,it5)),...
- 'color',(1-alpha)*[0.8500, 0.3250, 0.0980]+alpha*[0, 0.4470, 0.7410],...
- 'DisplayName',['t=',num2str(Ts5D(it5))]); hold on;
- end
- grid on; ylim([1e0,1e10]);
- xlabel('$p$');
- TITLE = ['$\sum_{kx,ky} |N_i^{p',num2str(Ji(ij)),'}|^2$']; title(TITLE);
-end
-save_figure
-end
-%%
-if 0
-%% Laguerre energy spectra
-% tf = Ts2D(end-3);
-skip = 1;
-fig = figure; FIGNAME = ['laguerre_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 1800, 600]);
-plt = @(x) squeeze(x);
-N10 = floor(Npi/10);
-Nrow = ceil((1+N10)/2);
-Ncol = ceil(Npi/Nrow/2);
-for ip = 1:2:Npi
- subplot(Nrow,Ncol,ip/2+1)
- for it5 = 1:skip:Ns5D
- alpha = it5*1.0/Ns5D;
- loglog(Ji,plt(epsilon_i_pj(ip,:,it5)),...
- 'color',(1-alpha)*[0.8500, 0.3250, 0.0980]+alpha*[0, 0.4470, 0.7410],...
- 'DisplayName',['t=',num2str(Ts5D(it5))]); hold on;
- grid on;
- xlabel('$j$'); ylim([1e-20,1e10]);
- TITLE = ['$\sum_{kx,ky} |N_i^{',num2str(Pi(ip)),'j}|^2$']; title(TITLE);
- end
-end
-save_figure
-end
-
-%%
-no_AA = (2:floor(2*Nkx/3));
-tKHI = 100;
-[~,itKHI] = min(abs(Ts2D-tKHI));
-after_KHI = (itKHI:Ns2D);
-if 0
-%% Phi frequency space time diagram at ky=ky(iky)
-ky_ = 0.0;
-[~,iky] = min(abs(ky-ky_));
-fig = figure; FIGNAME = ['phi_freq_diag_',PARAMS];set(gcf, 'Position', [100, 100, 500, 400]);
- [TY,TX] = meshgrid(Ts2D(after_KHI),kx(no_AA));
- pclr = pcolor(TX,TY,(squeeze(abs(PHI(no_AA,iky,(after_KHI)))))); set(pclr, 'edgecolor','none'); colorbar;
- caxis([0,10000]); colormap hot
- ylabel('$t c_s/R$'), xlabel('$0<k_r<2/3 k_r^{\max}$')
- legend(['$|\tilde\phi(k_z=',num2str(ky_),')|$'])
- title('Spectrogram of $\phi$')
-end
\ No newline at end of file
diff --git a/wk/transport_results_2_5.m b/wk/transport_results_2_5.m
deleted file mode 100644
index df40b342799186f9ad0761c49971f19057f883f0..0000000000000000000000000000000000000000
--- a/wk/transport_results_2_5.m
+++ /dev/null
@@ -1,156 +0,0 @@
-default_plots_options
-%% nuDGGK = 1.0
-figure; set(gcf, 'Position', [100, 100, 1200, 350])
-% 6,3 nuDGGK=1.0
-eta = [0.5,0.6,0.7,0.8];
-gamma = [32.6443 3.6895 0.3744 0.056];
-std_g = [6.1529 0.7986 0.0493 0.0134];
-subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=1.0$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(1,:)); hold on;
-subplot(122)
-
-% 10,5 nuDGGK = 1.0
-eta = [0.5 0.6,0.7,0.8];
-gamma = [32.6 3.5546 0.3917 0.0500];
-std_g = [7.7 0.5846 0.0486 0.0088];
-shear = [1.8505 0.60866 0.048249];
-std_s = [0.1599 0.00614 0.0061403];
-subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=1.0$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(2,:)); hold on;
-subplot(122)
-errorbar(eta(2:end),shear,std_s,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=1.0$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(2,:)); hold on;
-
-% 12,6 nuDGGK = 1.0
-eta = [0.6];
-gamma = [4.064];
-std_g = [0.7964];
-subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=12,06$, $\nu_{DGGK}=1.0$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(3,:)); hold on;
-subplot(122)
-
-
-% Mix len Ricci 2006
-subplot(121)
-semilogy([0.5 1.0],[10 1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$')
-subplot(122)
-grate = [0.2131 0.106 0.02021];
-semilogy([0.6 0.7 0.8],grate,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$V_E$''')
-
-%% nuDGGK = 0.5
-figure; set(gcf, 'Position', [100, 100, 1200, 350])
-% 6,3 nuDGGK=0.5
-eta = [0.5,0.6,0.7,0.8];
-gamma = [20.511 2.6292 0.2353 0.057];
-std_g = [3.67 1.2 0.055 0.004];
-shear = [nan 1.7417 0.57345 0.25155];
-std_s = [nan 0.35991 0.041 0.00913];
-subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.5$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(4,:)); hold on;
-subplot(122)
-errorbar(eta,shear,std_s,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.5$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(4,:)); hold on;
-
-% 10,5 nuDGGK=0.5
-eta = [0.6,0.7,0.8 0.9];
-gamma = [2.4949 0.23368 0.057792 0.023572];
-std_g = [0.8696 0.085267 0.0060463 0.0046137];
-shear = [1.7094 0.55278 0.2054 0.085678];
-std_s = [0.2428 0.068223 0.012734 0.012291];
-subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.5$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(5,:)); hold on;
-subplot(122)
-errorbar(eta,shear,std_s,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.5$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(5,:)); hold on;
-
-% Mix len Ricci 2006
-subplot(121)
-semilogy([0.5 1.0],[10 1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$')
-subplot(122)
-grate = [0.2194 0.129 0.05084 0.01346];
-semilogy([0.6 0.7 0.8 0.9],grate,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$V_E$''')
-
-%% nuDGGK = 0.1
-figure; set(gcf, 'Position', [100, 100, 1200, 350])
-% 6,3
-eta = [0.6 0.7 0.8];
-gamma = [0.24321 0.085 0.0367];
-std_g = [0.295 0.05 0.0023];
-gbmax = [1.0 0.21 0.0367];
-gbmin = [0.02 0.04 0.0367];
-% errorbar(eta,gamma,std_,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.1$',...
-% 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
-subplot(121)
-plot(eta,gamma,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
-plot(eta,gbmax,'-.','DisplayName',' ','Color',line_colors(6,:)); hold on;
-plot(eta,gbmin,'-.','DisplayName',' ','Color',line_colors(6,:)); hold on;
-
-% 10,5
-eta = [0.5 0.6 0.7 0.8 0.9];
-gamma = [12.2 0.19761 0.088 0.04253 0.026037];
-std_g = [4.7 0.21328 0.065 0.009 0.00118];
-gbmax = [12.2 0.8 0.25 0.04253 0.026037];
-gbmin = [12.2 0.02 0.03 0.04253 0.026037];
-shear = [0.65361 0.46548 0.30645 0.21123 ];
-std_s = [0.21288 0.10233 0.02886 0.0044664 ];
-sbmax = [1 0.7 0.30645 0.21123 ];
-sbmin = [0.4 0.35 0.30645 0.21123 ];
-% errorbar(eta,gamma,std_,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
-% 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
-subplot(121)
-plot(eta,gamma,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
-plot(eta,gbmax,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
-plot(eta,gbmin,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
-subplot(122)
-plot(eta(2:end),shear,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
-plot(eta(2:end),sbmax,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
-plot(eta(2:end),sbmin,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
-
-% Mix len Ricci 2006
-subplot(121)
-semilogy([0.5 1.0],[10 1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$')
-subplot(122)
-grate = [0.236 0.174 0.112 0.053];
-semilogy([0.6 0.7 0.8 0.9],grate,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$V_E$''')
-
-%% nuSGGK = 1.0
-figure
-% 6,3 nuDGGK=1.0
-eta = [0.5 0.6,0.7,0.8];
-gamma = [2.3 0.2215 0.0133 0.0032];
-std_g = [3.1 0.22 0.0019 0.0006];
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{SGGK}=1.0$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(8,:)); hold on;
-
-
-% 10,5 nuDGGK = 1.0
-eta = [0.5 0.6,0.7,0.8];
-gamma = [10 0.319 0.009 0.0026];
-std_g = [1.34 0.228 0.001 0.001];
-% errorbar(eta,gamma,err,'-.','DisplayName','$P,J=10,5$, $\nu_{SGGK}=1.0$',...
-% 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(1,:)); hold on;
-
-% Mix len Ricci 2006
-semilogy([0.5 1.0],[10 1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$')
diff --git a/wk/transport_results_2_6.m b/wk/transport_results_2_6.m
deleted file mode 100644
index b329b17403c51339512d9d3c488628eac2396157..0000000000000000000000000000000000000000
--- a/wk/transport_results_2_6.m
+++ /dev/null
@@ -1,146 +0,0 @@
-default_plots_options
-%% nuDGGK = 0.1
-figure; set(gcf, 'Position', [100, 100, 1200, 350])
-% 6,3 nuDGGK=0.1
-eta = [0.5,0.6,0.7];
-gamma = [116 24.396 3.684];
-std_g = [23 3 0.4];
-subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.1$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
-subplot(122)
-
-% 10,5 nuDGGK = 0.1
-eta = [0.5 0.6,0.7];
-gamma = [110 24.9 4.1];
-std_g = [23 4 0.38];
-shear = [4.6 4.0 1.35];
-std_s = [0.5 0.41 0.12];
-subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
-subplot(122)
-errorbar(eta,shear,std_s,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
-
-kmax = [0.576 0.5236 0.314];
-gmax = [0.305 0.20 0.06];
-subplot(121)
-% Mix len Ricci 2006
-semilogy([0.5 1.0],[10 1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
-% Bohm transport :
-btransp = gmax./kmax.^2;
-semilogy(.5:.1:.7,btransp','--','color',[0.4,0,0]+0.6, 'DisplayName','$\gamma/k^2$');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$')
-subplot(122)
-semilogy([0.5 0.6 0.7],gmax,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$V_E$''')
-
-%% nuDGGK = 0.01
-figure; set(gcf, 'Position', [100, 100, 1200, 350])
-% 6,3 nuDGGK=0.01
-eta = [0.5,0.6,0.7,0.8];
-gamma = [30.2 4.3 0.37 0.06];
-std_g = [4 0.7 0.05 0.008];
-shear = [5.0 2.0 0.6 0.16];
-std_s = [0.37 0.17 0.03 0.012];
-subplot(121)
-errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.01$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', 'b'); hold on;
-subplot(122)
-errorbar(eta,shear,std_s,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.01$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', 'b'); hold on;
-
-% % 10,5 nuDGGK=0.01
-% eta = [0.6,0.7,0.8 0.9];
-% gamma = [2.4949 0.23368 0.057792 0.023572];
-% std_g = [0.8696 0.085267 0.0060463 0.0046137];
-% shear = [1.7094 0.55278 0.2054 0.085678];
-% std_s = [0.2428 0.068223 0.012734 0.012291];
-% subplot(121)
-% errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.5$',...
-% 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(5,:)); hold on;
-% subplot(122)
-% errorbar(eta,shear,std_s,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.5$',...
-% 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(5,:)); hold on;
-
-kmax = [0.733 0.63 0.52 0.47];
-gmax = [0.32 0.22 0.11 0.027];
-subplot(121)
-% Mix len Ricci 2006
-semilogy([0.5 1.0],[10 1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
-% Bohm transport :
-btransp = gmax./kmax.^2;
-semilogy(.5:.1:.8,btransp','--','color',[0.4,0,0]+0.6, 'DisplayName','$\gamma/k^2$');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$')
-subplot(122)
-semilogy([0.5:0.1:0.8],gmax,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$V_E$''')
-
-%% nuDGGK = 0.001
-figure; set(gcf, 'Position', [100, 100, 1200, 350])
-% 6,3
-eta = [0.6 0.7 0.8 0.9];
-gamma = [0.26 0.088 0.042 0.0156];
-std_g = [0.4 0.06 0.003 0.0014];
-gbmax = [1.2 0.33 0.042 0.0156];
-gbmin = [0.015 0.035 0.003 0.0156];
-shear = [0.75 0.55 0.34 0.19];
-std_s = [0.33 0.12 0.012];
-sbmax = [1.4 0.796 0.34 0.19];
-sbmin = [0.4 0.4 0.34 0.19];
-
-subplot(121)
-plot(eta,gamma,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.001$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', 'k'); hold on;
-subplot(122)
-plot(eta,shear,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.001$',...
- 'LineWidth',2.0,'MarkerSize',5, 'Color', 'k'); hold on;
-
-% 10,5
-% eta = [0.5 0.6 0.7 0.8 0.9];
-% gamma = [12.2 0.19761 0.088 0.04253 0.026037];
-% std_g = [4.7 0.21328 0.065 0.009 0.00118];
-% gbmax = [12.2 0.8 0.25 0.04253 0.026037];
-% gbmin = [12.2 0.02 0.03 0.04253 0.026037];
-% shear = [0.65361 0.46548 0.30645 0.21123 ];
-% std_s = [0.21288 0.10233 0.02886 0.0044664 ];
-% sbmax = [1 0.7 0.30645 0.21123 ];
-% sbmin = [0.4 0.35 0.30645 0.21123 ];
-% errorbar(eta,gamma,std_,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
-% 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
-% subplot(121)
-% plot(eta,gamma,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
-% 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
-% plot(eta,gbmax,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
-% plot(eta,gbmin,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
-% subplot(122)
-% plot(eta(2:end),shear,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
-% 'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
-% plot(eta(2:end),sbmax,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
-% plot(eta(2:end),sbmin,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
-
-kmax = [0.27 0.68 0.52 0.31];
-gmax = [0.838 0.195 0.109 0.029];
-subplot(121)
-% Mix len Ricci 2006
-semilogy([0.5 1.0],[10 1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
-% Bohm transport :
-btransp = gmax./kmax.^2;
-semilogy(.6:.1:.9,btransp','--','color',[0.4,0,0]+0.6, 'DisplayName','$\gamma/k^2$');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$')
-subplot(122)
-semilogy([.6:.1:.9],gmax,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
-set(gca, 'YScale', 'log'); grid on; legend('show')
-xlabel('$L_n/R$'); ylabel('$V_E$''')
-subplot(121)
-plot(eta,gbmax,'-.','DisplayName',' ','Color','k'); hold on;
-plot(eta,gbmin,'-.','DisplayName',' ','Color','k'); hold on;
-subplot(122)
-plot(eta,sbmax,'-.','DisplayName',' ','Color','k'); hold on;
-plot(eta,sbmin,'-.','DisplayName',' ','Color','k'); hold on;
diff --git a/wk/untitled.m b/wk/untitled.m
deleted file mode 100644
index 73a439493a4e5880f160dbb13be1d294a151769c..0000000000000000000000000000000000000000
--- a/wk/untitled.m
+++ /dev/null
@@ -1,32 +0,0 @@
-N = 100;
-L = 20;
-
-A0 = 2;
-N0 = 4;
-K0 = 2*pi/L*N0;
-
-f = @(x_) A0 * sin(K0*x_);
-
-
-x = linspace(0,L, N);
-
-dk= 2*pi/L;
-k = dk*(-N/2:N/2-1);
-
-F = zeros(1,N);
-
-[~,ik0p] = min(abs(k-K0));
-[~,ik0m] = min(abs(k+K0));
-
-F(ik0p) = -A0/2 * 1i/dk;
-F(ik0m) = A0/2 * 1i/dk;
-
-
-f_ = ifft(fftshift(F));
-figure
-plot(x,f(x)); hold on;
-plot(x,N*dk*f_)
-
-%%
-figure
-plot(k,imag(F))
\ No newline at end of file