Script to compare simple simulations with MOLI

matlab/MOLI_time_solver.m

@@ -25,7 +25,7 @@ options.MOLI = 1;    % 1 -> Solve MOLI, 0 -> off
 options.solver.solver = 3;
 
 % MOLI Linear Fit Solver
-options.LinFitSolver = 2;
+options.LinFitSolver = 0;
 
 %% Main parameter scan
 
@@ -86,16 +86,16 @@ options.COSOlver.cmd   = 'mpirun -np 6 ./bin/CO 2 2 2';
 options.magnetic        = 1;           % 1-> Add toroidal magnetic gradient drift resonance effects
 
 % Physical Parameters
-params.tau              = 1.0;         % 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.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.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               = 1.0;	       % Major Radius / Background density gradient length
-params.RTe              = 0.0;         % Major Radius * normalized kperp / Background electron temperature gradient length
-params.RTi              = 0.0;         % Major Radius * normalized kperp / Background ion temperature gradient length
+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.
 
@@ -143,7 +143,6 @@ 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] ...
-params.RK4 = 1;
 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

wk/benchmark.m

+SIMID = 'benchmark'; % 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 = 20;
+GRID.jmaxe = 10;
+GRID.pmaxi = 20;
+GRID.jmaxi = 10;
+GRID.nkr   = 1;
+GRID.krmin = 0.;
+GRID.krmax = 0.;
+GRID.nkz   = 1;
+GRID.kzmin = 1.0;
+GRID.kzmax = 1.0;
+%% Model parameters
+MODEL.CO      = -1;  % Collision operator (-1 = Full Coulomb, 0 = Dougherty)
+MODEL.nu      = 0.0; % 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   = 0.0;        % Density
+MODEL.eta_T   = 1.0;        % Temperature
+MODEL.eta_B   = 0.0;        % 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_time_solver
+if params.RK4; MOLIsolvername = 'RK4';
+else;          MOLIsolvername = 'ode15i';
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Analysis and basic figures
+SAVEFIG = 1;
+filename = 'results_00.h5';
+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),')$'];
+%% Nipj
+
+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
+
+%% phi
+timephi  = h5read(filename,'/data/var2d/time');
+kr       = h5read(filename,'/data/var2d/phi/coordkr');
+kz       = h5read(filename,'/data/var2d/phi/coordkz');
+phiHeLaZ      = zeros(numel(timephi),numel(kr),numel(kz));
+for it = 1:numel(timephi)
+    tmp         = h5read(filename,['/data/var2d/phi/' num2str(it,'%06d')]);
+    phiHeLaZ(it,:,:) = tmp.real + 1i * tmp.imaginary;
+end
+
+timephiMOLI = results.time;
+phiMOLI  = zeros(size(timephiMOLI));
+for it = 1:numel(timephiMOLI)
+    phiMOLI(it) = get_phi(results.Napj(it,:),params,options);
+end
+
+fig = figure;
+%HeLaZ results
+semilogy(timephi,abs(phiHeLaZ),'-','DisplayName','HeLaZ RK4')
+hold on
+title(TITLE);
+%MOLI results
+semilogy(timephiMOLI,abs(phiMOLI),'--','DisplayName',['MOLI ',MOLIsolvername])
+grid on
+xlabel('$t$')
+ylabel('$|\phi|$')
+legend('show')
+if SAVEFIG; FIGNAME = 'phi'; save_figure; end;
+
+
+fig = figure;
+%HeLaZ results
+x1 = timeNi;
+y1 = abs(Nipj);
+semilogy(x1,y1,'-','DisplayName','HeLaZ RK4')
+hold on
+%MOLI results
+x2 = results.time;
+y2 = abs(results.Napj(:,1));
+semilogy(x2,y2,'--','DisplayName',['MOLI ',MOLIsolvername]);
+title(TITLE);
+grid on
+legend('show')
+xlabel('$t$')
+ylabel(['$|',moment,'|$'])
+if SAVEFIG; FIGNAME = 'Ni00'; save_figure; end;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file