diff --git a/.gitignore b/.gitignore index 500b87e..c6daf83 100644 --- a/.gitignore +++ b/.gitignore @@ -1,30 +1,31 @@ .vimrc *.swp *~ *.asv *.m.orig *.DS_Store *.fig *.pdf *.eps *.mov *.mp4 *.emf *.pptx *.jpg *.jpeg *.mat *.dat *.mod *.h5 *.o *.a srcinfo srcinfo.h logs/ results/ mod/ obj/ bin/ .gitignore wk/fort.90 +.directory diff --git a/wk/fort.90 b/wk/fort.90 index 73a6256..84051bf 100644 --- a/wk/fort.90 +++ b/wk/fort.90 @@ -1,52 +1,52 @@ &BASIC nrun=100000 dt=0.01 - tmax=200 ! time normalized to 1/omega_pe + tmax=50 ! time normalized to 1/omega_pe / &GRID - pmaxe =15 ! number of Hermite moments - jmaxe = 4 ! number of Hermite moments - pmaxi = 15 ! number of Hermite moments - jmaxi = 4 ! number of Hermite moments + pmaxe =81 ! number of Hermite moments + jmaxe = 20 ! number of Hermite moments + pmaxi = 81 ! number of Hermite moments + jmaxi = 20 ! number of Hermite moments nkr = 1 krmin = 0 krmax = 0 ! Normalized to cs0/Omega_i - nkz = 9 - kzmin = 0.1 + nkz = 1 + kzmin = 1 kzmax = 1 ! Normalized to cs0/Omega_i / &OUTPUT_PAR nsave_0d = 0 nsave_1d = 0 nsave_2d = 100 nsave_5d = 0 write_Ni00 = .true. write_moments = .true. write_phi = .true. write_doubleprecision = .true. resfile0 = 'results' / &MODEL_PAR ! Collisionality nu = 0.001 ! Normalized collision frequency normalized to plasma frequency tau_e = 1 ! T_e/T_e tau_i = 1 ! T_i/T_e temperature ratio sigma_e = 0.023338 ! sqrt(m_e/m_i) mass ratio sigma_i = 1 ! sqrt(m_i/m_i) q_e = -1 ! Electrons charge q_i = 1 ! Ions charge eta_n = 1 ! L_perp / L_n Density gradient eta_T = 0 ! L_perp / L_T Temperature gradient eta_B = 0.5 ! L_perp / L_B Magnetic gradient and curvature lambdaD = 0 ! Debye length / &INITIAL_CON ! Background value initback_moments=0.01 ! x 1e-3 ! Noise amplitude initnoise_moments=0 iseed=42 / &TIME_INTEGRATION_PAR numerical_scheme='RK4' / \ No newline at end of file diff --git a/wk/launcher.m b/wk/launcher.m index fd1c83d..9f031c5 100644 --- a/wk/launcher.m +++ b/wk/launcher.m @@ -1,64 +1,141 @@ -SIMID = 'Linear_Zpinch'; % Name of the simulations +SIMID = 'MOLI_Comparison'; % Name of the simulations addpath(genpath('../matlab')) % ... add %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 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 = 15; -GRID.jmaxe = 4; -GRID.pmaxi = 15; -GRID.jmaxi = 4; +GRID.pmaxe = 25; +GRID.jmaxe = 10; +GRID.pmaxi = 25; +GRID.jmaxi = 10; GRID.nkr = 1; GRID.krmin = 0.; GRID.krmax = 0.; -GRID.nkz = 9; +GRID.nkz = 1; GRID.kzmin = 0.1; -GRID.kzmax = 1.0; +GRID.kzmax = 0.1; %% Model parameters -MODEL.nu = 0.001; +MODEL.nu = 0.1; 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 = 200; +BASIC.nrun = 100000; +BASIC.dt = 0.01; +BASIC.tmax = 200; 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 solver +%% Run HeLaZ nproc = 1; EXEC = ' ../bin/helaz '; RUN = ['mpirun -np ' num2str(nproc)]; CMD = [RUN, EXEC, INPUT]; system(CMD); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Run MOLI +MOLI_time_solver + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Analysis and basic figures SAVEFIG = 1; -helaz_analysis +filename = 'results_00.h5'; +%% 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 + + +%% Error +nsave = OUTPUTS.nsave_2d; +if(numel(phiMOLI(1:nsave:end)) == numel(phiHeLaZ)) + errphi = abs(phiMOLI(1:nsave:end)-phiHeLaZ)./abs(phiMOLI(1:nsave:end)); + errNipj = abs(results.NapjRK4(1:nsave:end,1)-Nipj)./abs(results.NapjRK4(1:nsave:end,1)); + figure + plot(timephi,errphi*100,'-','DisplayName','$\epsilon(\phi)$') + hold on; + plot(timephi,errNipj*100,'--','DisplayName','$\epsilon(N_i^{00})$') + title(TITLE); + xlabel('$t$'); + ylabel('$\epsilon$ [\%]') + grid on + legend('show') +else + figure + %HeLaZ results + plot(timephi,abs(phiHeLaZ),'-','DisplayName','HeLaZ RK4') + title(TITLE); + hold on + %MOLI results + plot(timephiMOLI,abs(phiMOLI),'--','DisplayName','MOLI RK4') + grid on + xlabel('$t$') + ylabel('$|\phi|$') + legend('show') + + figure + %HeLaZ results + x1 = timeNi; + y1 = abs(Nipj); + plot(x1,y1,'-','DisplayName','HeLaZ RK4') + hold on + %MOLI results + x2 = results.timeRK4; + y2 = abs(results.NapjRK4(:,1)); + plot(x2(1:100:end),y2(1:100:end),'--','DisplayName','MOLI RK4'); + title(TITLE); + grid on + legend('show') + xlabel('$t$') + ylabel(['$|',moment,'|$']) +end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \ No newline at end of file