From 43291f7e97b12c99f70dcb3bdc1187cfe06ef328 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Fri, 3 Jul 2020 18:23:43 +0200
Subject: [PATCH] tests for kperp scan
---
matlab/MOLI_kperp_scan.m | 164 ++++++++++++++++++++++++++++++++++++++
wk/benchmark_kperp_scan.m | 134 +++++++++++++++++++++++++++++++
2 files changed, 298 insertions(+)
create mode 100644 matlab/MOLI_kperp_scan.m
create mode 100644 wk/benchmark_kperp_scan.m
diff --git a/matlab/MOLI_kperp_scan.m b/matlab/MOLI_kperp_scan.m
new file mode 100644
index 00000000..8ed91498
--- /dev/null
+++ b/matlab/MOLI_kperp_scan.m
@@ -0,0 +1,164 @@
+%% 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 = 2;
+
+% 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 = GRID.kzmin; % normalized perpendicular wave number to the soundLarmor radius. Note: If ions ==0 (e.g. EPW), kperp --> b
+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.initback_moments; % 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/wk/benchmark_kperp_scan.m b/wk/benchmark_kperp_scan.m
new file mode 100644
index 00000000..bdec0eb9
--- /dev/null
+++ b/wk/benchmark_kperp_scan.m
@@ -0,0 +1,134 @@
+SIMID = 'benchmark_kperp_scan'; % Name of the simulations
+addpath(genpath('../matlab')) % ... add
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% outputs options
+OUTPUTS.nsave_0d = 0;
+OUTPUTS.nsave_1d = 0;
+OUTPUTS.nsave_2d = 1;
+OUTPUTS.nsave_5d = 0;
+OUTPUTS.write_Ni00 = '.true.';
+OUTPUTS.write_moments = '.true.';
+OUTPUTS.write_phi = '.true.';
+OUTPUTS.write_doubleprecision = '.true.';
+OUTPUTS.resfile0 = '''results''';
+%% Grid parameters
+GRID.pmaxe = 15;
+GRID.jmaxe = 6;
+GRID.pmaxi = 15;
+GRID.jmaxi = 6;
+GRID.nkr = 1;
+GRID.krmin = 0.;
+GRID.krmax = 0.;
+GRID.nkz = 20;
+GRID.kzmin = 0.1;
+GRID.kzmax = 1.5;
+%% Model parameters
+MODEL.CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+MODEL.nu = 0.1; % collisionality nu*L_perp/Cs0
+% temperature ratio T_a/T_e
+MODEL.tau_e = 1.0;
+MODEL.tau_i = 1.0;
+% mass ratio sqrt(m_a/m_i)
+MODEL.sigma_e = 0.0233380;
+MODEL.sigma_i = 1.0;
+% charge q_a/e
+MODEL.q_e =-1.0;
+MODEL.q_i = 1.0;
+% gradients L_perp/L_x
+MODEL.eta_n = 1.0; % Density
+MODEL.eta_T = 0.0; % Temperature
+MODEL.eta_B = 0.5; % Magnetic
+% Debye length
+MODEL.lambdaD = 0.0;
+%% Time integration and intialization parameters
+TIME_INTEGRATION.numerical_scheme = '''RK4''';
+BASIC.nrun = 100000;
+BASIC.dt = 0.05;
+BASIC.tmax = 100.0;
+INITIAL.initback_moments = 0.01;
+INITIAL.initnoise_moments = 0.;
+INITIAL.iseed = 42;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Write input file
+INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Run HeLaZ
+nproc = 1;
+EXEC = ' ../bin/helaz ';
+RUN = ['mpirun -np ' num2str(nproc)];
+CMD = [RUN, EXEC, INPUT];
+system(CMD);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Run MOLI with same input parameters
+params.RK4 = 1;
+run ../matlab/MOLI_kperp_scan
+if params.RK4; MOLIsolvername = 'RK4';
+else; MOLIsolvername = 'ode15i';
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Analysis and basic figures
+SAVEFIG = 1;
+filename = 'results_00.h5';
+
+if MODEL.CO == -1; CONAME = 'FC';
+elseif MODEL.CO == -2; CONAME = 'DC';
+elseif MODEL.CO == 0; CONAME = 'LB'; end;
+
+TITLE = [];
+TITLE = [TITLE,'$\eta_n=',num2str(MODEL.eta_n),'$, '];
+TITLE = [TITLE,'$\eta_B=',num2str(MODEL.eta_B),'$, '];
+TITLE = [TITLE,'$\eta_T=',num2str(MODEL.eta_T),'$, '];
+TITLE = [TITLE, '$\nu=',num2str(MODEL.nu),'$, '];
+TITLE = [TITLE, '$(P,J)=(',num2str(GRID.pmaxe),',',num2str(GRID.jmaxe),')$'];
+
+%% Growth rate analysis
+gammas = zeros(numel(kr),numel(kz));
+shifts = zeros(numel(kr),numel(kz));
+
+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
+
+% Linear fit of log(Napj)
+x1 = timeNi;
+itmin = ceil(0.5 * numel(timeNi)); %Take the second half of the time evolution
+ikr = 1;
+
+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
+
+%% Plot
+fig = figure;
+
+%HeLaZ results
+X = kz*sqrt(1+MODEL.tau_i);
+Y = gammas(1,:);
+plot(X,Y,'-','DisplayName','HeLaZ')
+hold on
+
+%MOLI results
+X = kz*sqrt(1+MODEL.tau_i);
+Y = results.kperp.Maxgammas;
+plot(X,Y,'--','DisplayName','MOLI');
+title(TITLE);
+grid on
+legend('show')
+xlabel('$k_\perp * (1+\tau)^{1/2}$')
+ylabel('$\gamma L_\perp/c_{s} $')
+if SAVEFIG; FIGNAME = 'g_vs_kz'; save_figure; end;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
--
GitLab