script to compare linear results with MOLI (RK4 needed in MOLI)

wk/MOLI_comparison.m

+SIMID = 'MOLI_Comparison'; % Name of the simulations
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% outputs options
+OUTPUTS.nsave_0d = 0;
+OUTPUTS.nsave_1d = 0;
+OUTPUTS.nsave_2d = 100;
+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 = 81;
+GRID.jmaxe = 20;
+GRID.pmaxi = 81;
+GRID.jmaxi = 20;
+GRID.nkr   = 1;
+GRID.krmin = 0.;
+GRID.krmax = 0.;
+GRID.nkz   = 1;
+GRID.kzmin = 1.0;
+GRID.kzmax = 1.0;
+%% Model parameters
+MODEL.nu      = 0.001;
+MODEL.tau_e   = 1.0;
+MODEL.tau_i   = 1.0;
+MODEL.sigma_e = 0.0233380;
+MODEL.sigma_i = 1.0;
+MODEL.q_e     =-1.0;
+MODEL.q_i     = 1.0;
+MODEL.eta_n   = 1.0;
+MODEL.eta_T   = 0.0;
+MODEL.eta_B   = 0.5;
+MODEL.lambdaD = 0.0;
+%% Time integration and intialization parameters
+TIME_INTEGRATION.numerical_scheme  = '''RK4''';
+BASIC.nrun                = 100000;
+BASIC.dt                  = 0.01;
+BASIC.tmax                = 50.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
+run ../matlab/MOLI_time_solver
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Analysis and basic figures
+SAVEFIG = 1;
+filename = 'results_00.h5';
+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, '$(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.timeRK4;
+phiMOLI  = zeros(size(timephiMOLI));
+for it = 1:numel(timephiMOLI)
+    phiMOLI(it) = get_phi(results.NapjRK4(it,:),params,options);
+end
+
+fig = figure;
+%HeLaZ results
+semilogy(timephi,abs(phiHeLaZ),'-','DisplayName','HeLaZ RK4')
+title(TITLE);
+hold on
+%MOLI results
+semilogy(timephiMOLI,abs(phiMOLI),'--','DisplayName','MOLI RK4')
+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.timeRK4;
+y2 = abs(results.NapjRK4(:,1));
+semilogy(x2(1:100:end),y2(1:100:end),'--','DisplayName','MOLI RK4');
+title(TITLE);
+grid on
+legend('show')
+xlabel('$t$')
+ylabel(['$|',moment,'|$'])
+if SAVEFIG; FIGNAME = 'Ni00'; save_figure; end;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file