diff --git a/Makefile b/Makefile index d776c7d..21a929d 100644 --- a/Makefile +++ b/Makefile @@ -1,156 +1,159 @@ include local/dirs.inc include local/make.inc EXEC = $(BINDIR)/helaz F90 = mpiifort # Add Multiple-Precision Library EXTLIBS += -L$(FMDIR)/lib EXTINC += -I$(FMDIR)/mod EXTLIBS += -L$(FFTWDIR)/lib EXTINC += -I$(FFTWDIR)/include all: dirs src/srcinfo.h $(EXEC) install: dirs src/srcinfo.h $(EXEC) mvmod run: all (cd wk; $(EXEC);) dirs: mkdir -p $(BINDIR) mkdir -p $(OBJDIR) mkdir -p $(MODDIR) src/srcinfo.h: ( cd src/srcinfo; $(MAKE);) clean: cleanobj cleanmod @rm -f src/srcinfo.h @rm -f src/srcinfo/srcinfo.h cleanobj: @rm -f $(OBJDIR)/*o cleanmod: @rm -f $(MODDIR)/*mod @rm -f *.mod cleanbin: @rm -f $(EXEC) mvmod: mv *.mod mod/. $(OBJDIR)/diagnose.o : src/srcinfo.h FOBJ=$(OBJDIR)/advance_field.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o $(OBJDIR)/basic_mod.o \ $(OBJDIR)/coeff_mod.o $(OBJDIR)/closure_mod.o $(OBJDIR)/collision_mod.o \ $(OBJDIR)/compute_Sapj.o $(OBJDIR)/control.o $(OBJDIR)/fourier_mod.o \ $(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o $(OBJDIR)/fields_mod.o \ $(OBJDIR)/inital.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/main.o $(OBJDIR)/memory.o \ $(OBJDIR)/model_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o \ -$(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/readinputs.o \ -$(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/restarts_mod.o $(OBJDIR)/stepon.o \ -$(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o +$(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/processing_mod.o \ +$(OBJDIR)/readinputs.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/restarts_mod.o \ +$(OBJDIR)/stepon.o $(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(EXEC): $(FOBJ) $(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@ $(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/advance_field.F90 -o $@ $(OBJDIR)/array_mod.o : src/array_mod.F90 $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/array_mod.F90 -o $@ $(OBJDIR)/auxval.o : src/auxval.F90 $(OBJDIR)/fourier_mod.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o $(OBJDIR)/grid_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/auxval.F90 -o $@ $(OBJDIR)/basic_mod.o : src/basic_mod.F90 $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/basic_mod.F90 -o $@ $(OBJDIR)/coeff_mod.o : src/coeff_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/coeff_mod.F90 -o $@ $(OBJDIR)/closure_mod.o : src/closure_mod.F90 $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/closure_mod.F90 -o $@ $(OBJDIR)/collision_mod.o : src/collision_mod.F90 $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/collision_mod.F90 -o $@ $(OBJDIR)/compute_Sapj.o : src/compute_Sapj.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/compute_Sapj.F90 -o $@ $(OBJDIR)/control.o : src/control.F90 $(OBJDIR)/auxval.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/readinputs.o $(OBJDIR)/tesend.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/control.F90 -o $@ - $(OBJDIR)/fourier_mod.o : src/fourier_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o - $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_mod.F90 -o $@ - - $(OBJDIR)/diagnose.o : src/diagnose.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o + $(OBJDIR)/diagnose.o : src/diagnose.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/diagnose.F90 -o $@ $(OBJDIR)/diagnostics_par_mod.o : src/diagnostics_par_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/diagnostics_par_mod.F90 -o $@ $(OBJDIR)/endrun.o : src/endrun.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/endrun.F90 -o $@ $(OBJDIR)/fields_mod.o : src/fields_mod.F90 $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fields_mod.F90 -o $@ + $(OBJDIR)/fourier_mod.o : src/fourier_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o + $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_mod.F90 -o $@ + $(OBJDIR)/ghosts_mod.o : src/ghosts_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/ppinit.o $(OBJDIR)/time_integration_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ghosts_mod.F90 -o $@ $(OBJDIR)/grid_mod.o : src/grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/grid_mod.F90 -o $@ $(OBJDIR)/inital.o : src/inital.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/restarts_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/inital.F90 -o $@ $(OBJDIR)/initial_par_mod.o : src/initial_par_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/initial_par_mod.F90 -o $@ $(OBJDIR)/main.o : src/main.F90 $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/main.F90 -o $@ $(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/grid_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/memory.F90 -o $@ $(OBJDIR)/model_mod.o : src/model_mod.F90 $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/model_mod.F90 -o $@ $(OBJDIR)/moments_eq_rhs.o : src/moments_eq_rhs.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs.F90 -o $@ $(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/poisson.F90 -o $@ $(OBJDIR)/ppexit.o : src/ppexit.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppexit.F90 -o $@ $(OBJDIR)/ppinit.o : src/ppinit.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppinit.F90 -o $@ $(OBJDIR)/prec_const_mod.o : src/prec_const_mod.F90 $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/prec_const_mod.F90 -o $@ + $(OBJDIR)/processing_mod.o : src/processing_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/basic_mod.o + $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/processing_mod.F90 -o $@ + $(OBJDIR)/readinputs.o : src/readinputs.F90 $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/readinputs.F90 -o $@ $(OBJDIR)/restarts_mod.o : src/restarts_mod.F90 $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/restarts_mod.F90 -o $@ $(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/stepon.F90 -o $@ $(OBJDIR)/tesend.o : src/tesend.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/tesend.F90 -o $@ $(OBJDIR)/time_integration_mod.o : src/time_integration_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/time_integration_mod.F90 -o $@ $(OBJDIR)/utility_mod.o : src/utility_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@ diff --git a/matlab/compile_results.m b/matlab/compile_results.m index 74bd637..44065bb 100644 --- a/matlab/compile_results.m +++ b/matlab/compile_results.m @@ -1,64 +1,80 @@ CONTINUE = 1; JOBNUM = 0; JOBFOUND = 0; Nipj_ = []; Nepj_ = []; Ni00_ = []; Ne00_ = []; PHI_ = []; Ts2D_ = []; Ts5D_ = []; Sipj_ = []; Sepj_ = []; Pe_old = 1e9; Pi_old = Pe_old; Je_old = Pe_old; Ji_old = Pe_old; while(CONTINUE) filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM); if exist(filename, 'file') == 2 % Load results of simulation #JOBNUM load_results % Check polynomials degrees - sz = size(Nepj); Pe_new= sz(1); Je_new= sz(2); - sz = size(Nipj); Pi_new= sz(1); Ji_new= sz(2); + Pe_new= numel(Pe); Je_new= numel(Je); + Pi_new= numel(Pi); Ji_new= numel(Ji); % If a degree is larger than previous job, put them in a larger array if (sum([Pe_new, Je_new, Pi_new, Ji_new]>[Pe_old, Je_old, Pi_old, Ji_old]) >= 1) - tmp = Nipj_; sz = size(tmp); - Nipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')'); - Nipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp; - tmp = Nepj_; sz = size(tmp); - Nepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')'); - Nepj_(1:Pe_old,1:Je_old,:,:,:) = tmp; - tmp = Sipj_; sz = size(tmp); - Sipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')'); - Sipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp; - tmp = Sepj_; sz = size(tmp); - Sepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')'); - Sepj_(1:Pe_old,1:Je_old,:,:,:) = tmp; + if W_NAPJ + tmp = Nipj_; sz = size(tmp); + Nipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')'); + Nipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp; + tmp = Nepj_; sz = size(tmp); + Nepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')'); + Nepj_(1:Pe_old,1:Je_old,:,:,:) = tmp; + end + if W_SAPJ + tmp = Sipj_; sz = size(tmp); + Sipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')'); + Sipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp; + tmp = Sepj_; sz = size(tmp); + Sepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')'); + Sepj_(1:Pe_old,1:Je_old,:,:,:) = tmp; + end end - - Nipj_ = cat(5,Nipj_,Nipj); - Nepj_ = cat(5,Nepj_,Nepj); - Ni00_ = cat(3,Ni00_,Ni00); - Ne00_ = cat(3,Ne00_,Ne00); - PHI_ = cat(3,PHI_,PHI); - Ts2D_ = cat(1,Ts2D_,Ts2D); - Ts5D_ = cat(1,Ts5D_,Ts5D); - - Sipj_ = cat(5,Sipj_,Sipj); - Sepj_ = cat(5,Sepj_,Sepj); + + if W_PHI || W_NA00 + Ts2D_ = cat(1,Ts2D_,Ts2D); + end + if W_PHI + PHI_ = cat(3,PHI_,PHI); + end + if W_NA00 + Ni00_ = cat(3,Ni00_,Ni00); + Ne00_ = cat(3,Ne00_,Ne00); + end + + if W_NAPJ || W_SAPJ + Ts5D_ = cat(1,Ts5D_,Ts5D); + end + if W_NAPJ + Nipj_ = cat(5,Nipj_,Nipj); + Nepj_ = cat(5,Nepj_,Nepj); + end + if W_SAPJ + Sipj_ = cat(5,Sipj_,Sipj); + Sepj_ = cat(5,Sepj_,Sepj); + end JOBFOUND = JOBFOUND + 1; LASTJOB = JOBNUM; elseif (JOBNUM > 20) CONTINUE = 0; disp(['found ',num2str(JOBFOUND),' results']); end JOBNUM = JOBNUM + 1; Pe_old = Pe_new; Je_old = Je_new; Pi_old = Pi_new; Ji_old = Ji_new; end Nipj = Nipj_; Nepj = Nepj_; Ts5D = Ts5D_; Ni00 = Ni00_; Ne00 = Ne00_; PHI = PHI_; Ts2D = Ts2D_; clear Nipj_ Nepj_ Ni00_ Ne00_ PHI_ Ts2D_ Ts5D_ Sipj = Sipj_; Sepj = Sepj_; clear Sipj_ Sepj_ JOBNUM = LASTJOB filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM); \ No newline at end of file diff --git a/matlab/load_0D_data.m b/matlab/load_0D_data.m new file mode 100644 index 0000000..7c2182f --- /dev/null +++ b/matlab/load_0D_data.m @@ -0,0 +1,9 @@ +function [ data, time, dt ] = load_0D_data( filename, variablename ) +%LOAD_0D_DATA load a 0D variable stored in a hdf5 result file from HeLaZ + time = h5read(filename,'/data/var0d/time'); + dt = h5readatt(filename,'/data/input','dt'); + cstart= h5readatt(filename,'/data/input','start_iframe0d'); + data = h5read(filename,['/data/var0d/',variablename]); + +end + diff --git a/matlab/load_2D_data.m b/matlab/load_2D_data.m index f03bbaa..47ef91a 100644 --- a/matlab/load_2D_data.m +++ b/matlab/load_2D_data.m @@ -1,19 +1,16 @@ -function [ data, kr, kz, time, dt ] = load_2D_data( filename, variablename ) +function [ data, time, dt ] = load_2D_data( filename, variablename ) %LOAD_2D_DATA load a 2D variable stored in a hdf5 result file from HeLaZ time = h5read(filename,'/data/var2d/time'); kr = h5read(filename,'/data/grid/coordkr'); kz = h5read(filename,'/data/grid/coordkz'); -% kr = h5read(filename,['/data/var2d/',variablename,'/coordkr']); -% kz = h5read(filename,['/data/var2d/',variablename,'/coordkz']); - dt = h5readatt(filename,'/data/input','dt'); cstart= h5readatt(filename,'/data/input','start_iframe2d'); data = zeros(numel(kr),numel(kz),numel(time)); for it = 1:numel(time) tmp = h5read(filename,['/data/var2d/',variablename,'/', num2str(cstart+it,'%06d')]); data(:,:,it) = tmp.real + 1i * tmp.imaginary; end end diff --git a/matlab/load_5D_data.m b/matlab/load_5D_data.m index 87edc9b..e36d767 100644 --- a/matlab/load_5D_data.m +++ b/matlab/load_5D_data.m @@ -1,22 +1,23 @@ -function [ data, p, j, kr, kz, time, dt ] = load_5D_data( filename, variablename ) +function [ data, time, dt ] = load_5D_data( filename, variablename ) %LOAD_5D_DATA load a 5D variable stored in a hdf5 result file from HeLaZ time = h5read(filename,'/data/var5d/time'); - p = h5read(filename,'/data/grid/coordp'); - j = h5read(filename,'/data/grid/coordj'); + if strcmp(variablename,'moments_e') || strcmp(variablename,'Sepj') + p = h5read(filename,'/data/grid/coordp_e'); + j = h5read(filename,'/data/grid/coordj_e'); + else + p = h5read(filename,'/data/grid/coordp_i'); + j = h5read(filename,'/data/grid/coordj_i'); + end kr = h5read(filename,'/data/grid/coordkr'); kz = h5read(filename,'/data/grid/coordkz'); -% p = h5read(filename,['/data/var5d/',variablename,'/coordp']); -% j = h5read(filename,['/data/var5d/',variablename,'/coordj']); -% kr = h5read(filename,['/data/var5d/',variablename,'/coordkr']); -% kz = h5read(filename,['/data/var5d/',variablename,'/coordkz']); dt = h5readatt(filename,'/data/input','dt'); cstart= h5readatt(filename,'/data/input','start_iframe5d'); data = zeros(numel(p),numel(j),numel(kr),numel(kz),numel(time)); for it = 1:numel(time) tmp = h5read(filename,['/data/var5d/', variablename,'/', num2str(cstart+it,'%06d')]); data(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary; end end \ No newline at end of file diff --git a/matlab/load_grid_data.m b/matlab/load_grid_data.m new file mode 100644 index 0000000..4318157 --- /dev/null +++ b/matlab/load_grid_data.m @@ -0,0 +1,9 @@ +function [ pe, je, pi, ji, kr, kz ] = load_grid_data( filename ) +%LOAD_GRID_DATA stored in a hdf5 result file from HeLaZ + pe = h5read(filename,'/data/grid/coordp_e'); + je = h5read(filename,'/data/grid/coordj_e'); + pi = h5read(filename,'/data/grid/coordp_i'); + ji = h5read(filename,'/data/grid/coordj_i'); + kr = h5read(filename,'/data/grid/coordkr'); + kz = h5read(filename,'/data/grid/coordkz'); +end \ No newline at end of file diff --git a/matlab/load_params.m b/matlab/load_params.m index 0e47f49..c176b70 100644 --- a/matlab/load_params.m +++ b/matlab/load_params.m @@ -1,45 +1,51 @@ CO = h5readatt(filename,'/data/input','CO'); ETAB = h5readatt(filename,'/data/input','eta_B'); ETAN = h5readatt(filename,'/data/input','eta_n'); ETAT = h5readatt(filename,'/data/input','eta_T'); PMAXI = h5readatt(filename,'/data/input','pmaxi'); JMAXI = h5readatt(filename,'/data/input','jmaxi'); PMAXE = h5readatt(filename,'/data/input','pmaxe'); JMAXE = h5readatt(filename,'/data/input','jmaxe'); NON_LIN = h5readatt(filename,'/data/input','NON_LIN'); NU = h5readatt(filename,'/data/input','nu')/0.532; NR = h5readatt(filename,'/data/input','nr'); NZ = h5readatt(filename,'/data/input','nz'); L = h5readatt(filename,'/data/input','Lr'); CLOS = h5readatt(filename,'/data/input','CLOS'); MU = str2num(filename(end-18:end-14)); +W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y'; +W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y'; +W_NA00 = h5readatt(filename,'/data/input','write_Na00') == 'y'; +W_NAPJ = h5readatt(filename,'/data/input','write_Napj') == 'y'; +W_SAPJ = h5readatt(filename,'/data/input','write_Sapj') == 'y'; + if NON_LIN == 'y' NON_LIN = 1; else NON_LIN = 0; end if (CO == 0); CONAME = 'LB'; elseif(CO == -1); CONAME = 'FC'; elseif(CO == -2); CONAME = 'DG'; elseif(CO == -3); CONAME = 'DGGK'; end if (CLOS == 0); CLOSNAME = 'Trunc.'; elseif(CLOS == 1); CLOSNAME = 'Clos. 1'; elseif(CLOS == 2); CLOSNAME = 'Clos. 2'; end if (PMAXE == PMAXI) && (JMAXE == JMAXI) degngrad = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE)]; else degngrad = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),... '_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)]; end degngrad = [degngrad,'_eta_',num2str(ETAB/ETAN),'_nu_%0.0e_',... CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e']; degngrad = sprintf(degngrad,[NU,MU]); if ~NON_LIN; degngrad = ['lin_',degngrad]; end resolution = [num2str(NR),'x',num2str(NZ/2),'_']; gridname = ['L_',num2str(L),'_']; PARAMS = [resolution,gridname,degngrad]; % BASIC.RESDIR = [SIMDIR,PARAMS,'/']; diff --git a/matlab/load_results.m b/matlab/load_results.m index 1cb6be5..fc41f9c 100644 --- a/matlab/load_results.m +++ b/matlab/load_results.m @@ -1,15 +1,33 @@ %% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM); disp(['Loading ',filename]) % Loading from output file CPUTIME = h5readatt(filename,'/data/input','cpu_time'); DT_SIM = h5readatt(filename,'/data/input','dt'); -[Nipj, Pi, Ji, kr, kz, Ts5D, dt5D] = load_5D_data(filename, 'moments_i'); -[Nepj, Pe, Je ] = load_5D_data(filename, 'moments_e'); -[Ni00, kr, kz, Ts2D, dt2D] = load_2D_data(filename, 'Ni00'); - Ne00 = load_2D_data(filename, 'Ne00'); -PHI = load_2D_data(filename, 'phi'); +[Pe, Je, Pi, Ji, kr, kz] = load_grid_data(filename); -% Sipj = load_5D_data(filename, 'Sipj'); -% Sepj = load_5D_data(filename, 'Sepj'); +if W_GAMMA + [ GGAMMA_RI, Ts0D, dt0D] = load_0D_data(filename, 'gflux_ri'); + PGAMMA_RI = load_0D_data(filename, 'pflux_ri'); +end + +if W_PHI + [ PHI, Ts2D, dt2D] = load_2D_data(filename, 'phi'); +end + +if W_NA00 + [Ni00, Ts2D, dt2D] = load_2D_data(filename, 'Ni00'); + Ne00 = load_2D_data(filename, 'Ne00'); +end + + +if W_NAPJ + [Nipj, Ts5D, dt5D] = load_5D_data(filename, 'moments_i'); + [Nepj ] = load_5D_data(filename, 'moments_e'); +end + +if W_SAPJ + [Sipj, Ts5D, dt5D] = load_5D_data(filename, 'Sipj'); + Sepj = load_5D_data(filename, 'Sepj'); +end diff --git a/matlab/setup.m b/matlab/setup.m index 6f42c2b..a0da7b8 100644 --- a/matlab/setup.m +++ b/matlab/setup.m @@ -1,122 +1,127 @@ %% ________________________________________________________________________ SIMDIR = ['../results/',SIMID,'/']; % Grid parameters GRID.pmaxe = PMAXE; % Electron Hermite moments GRID.jmaxe = JMAXE; % Electron Laguerre moments GRID.pmaxi = PMAXI; % Ion Hermite moments GRID.jmaxi = JMAXI; % Ion Laguerre moments GRID.Nr = N; % r grid resolution GRID.Lr = L; % r length GRID.Nz = N * (1-KREQ0) + KREQ0; % z '' GRID.Lz = L * (1-KREQ0); % z '' GRID.kpar = KPAR; % Model parameters MODEL.CO = CO; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) MODEL.CLOS = CLOS; if NON_LIN; MODEL.NON_LIN = '.true.'; else; MODEL.NON_LIN = '.false.';end; MODEL.mu = MU; MODEL.mu_p = MU_P; MODEL.mu_j = MU_J; MODEL.nu = NU; % hyper diffusive coefficient nu for HW % temperature ratio T_a/T_e MODEL.tau_e = TAU; MODEL.tau_i = TAU; % 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; if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end; % gradients L_perp/L_x MODEL.eta_n = ETAN; % source term kappa for HW MODEL.eta_T = ETAT; % Temperature MODEL.eta_B = ETAB; % Magnetic MODEL.lambdaD = LAMBDAD; % if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KR0KH),'_A_',num2str(A0KH)]; end; % Time integration and intialization parameters TIME_INTEGRATION.numerical_scheme = '''RK4'''; if INIT_PHI; INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end; INITIAL.initback_moments = 0.0e-5; INITIAL.initnoise_moments = NOISE0; INITIAL.iseed = 42; fcmat_pmaxe = 25; fcmat_jmaxe = 18; fcmat_pmaxi = 25; fcmat_jmaxi = 18; INITIAL.selfmat_file = ... ['''../../../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_',num2str(fcmat_pmaxe),... '_Jmaxe_',num2str(fcmat_jmaxe),'_Pmaxi_',num2str(fcmat_pmaxi),'_Jmaxi_',... num2str(fcmat_jmaxi),'_pamaxx_10.h5''']; INITIAL.eimat_file = ... ['''../../../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_',num2str(fcmat_pmaxe),... '_Jmaxe_',num2str(fcmat_jmaxe),'_Pmaxi_',num2str(fcmat_pmaxi),'_Jmaxi_',... num2str(fcmat_jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5''']; INITIAL.iemat_file = ... ['''../../../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_',num2str(fcmat_pmaxe),... '_Jmaxe_',num2str(fcmat_jmaxe),'_Pmaxi_',num2str(fcmat_pmaxi),'_Jmaxi_',... num2str(fcmat_jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5''']; % Naming and creating input file if (CO == 0); CONAME = 'LB'; elseif(CO == -1); CONAME = 'FC'; elseif(CO == -2); CONAME = 'DG'; elseif(CO == -3); CONAME = 'DGGK'; end if (CLOS == 0); CLOSNAME = 'Trunc.'; elseif(CLOS == 1); CLOSNAME = 'Clos. 1'; elseif(CLOS == 2); CLOSNAME = 'Clos. 2'; end if (PMAXE == PMAXI) && (JMAXE == JMAXI) degngrad = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE)]; else degngrad = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),... '_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)]; end degngrad = [degngrad,'_eta_',num2str(ETAB/ETAN),'_nu_%0.0e_',... CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e']; degngrad = sprintf(degngrad,[NU,MU]); if ~NON_LIN; degngrad = ['lin_',degngrad]; end resolution = [num2str(GRID.Nr),'x',num2str(GRID.Nz/2),'_']; gridname = ['L_',num2str(L),'_']; PARAMS = [resolution,gridname,degngrad]; BASIC.RESDIR = [SIMDIR,PARAMS,'/']; BASIC.PARAMS = PARAMS; BASIC.SIMID = SIMID; BASIC.nrun = 1e8; BASIC.dt = DT; BASIC.tmax = TMAX; %time normalized to 1/omega_pe BASIC.maxruntime = str2num(CLUSTER.TIME(1:2))*3600 ... + str2num(CLUSTER.TIME(4:5))*60 ... + str2num(CLUSTER.TIME(7:8)); % Outputs parameters if RESTART; BASIC.RESTART = '.true.'; else; BASIC.RESTART = '.false.';end; OUTPUTS.nsave_0d = floor(1.0/SPS0D/DT); OUTPUTS.nsave_1d = -1; OUTPUTS.nsave_2d = floor(1.0/SPS2D/DT); OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT); OUTPUTS.nsave_cp = floor(1.0/SPSCP/DT); -OUTPUTS.write_doubleprecision = '.false.'; -OUTPUTS.resfile0 = '''outputs'''; -OUTPUTS.rstfile0 = '''checkpoint'''; -OUTPUTS.job2load = JOB2LOAD; +if W_DOUBLE; OUTPUTS.write_doubleprecision = '.true.'; else; OUTPUTS.write_doubleprecision = '.false.';end; +if W_GAMMA; OUTPUTS.write_gamma = '.true.'; else; OUTPUTS.write_gamma = '.false.';end; +if W_PHI; OUTPUTS.write_phi = '.true.'; else; OUTPUTS.write_phi = '.false.';end; +if W_NA00; OUTPUTS.write_Na00 = '.true.'; else; OUTPUTS.write_Na00 = '.false.';end; +if W_NAPJ; OUTPUTS.write_Napj = '.true.'; else; OUTPUTS.write_Napj = '.false.';end; +if W_SAPJ; OUTPUTS.write_Sapj = '.true.'; else; OUTPUTS.write_Sapj = '.false.';end; +OUTPUTS.resfile0 = '''outputs'''; +OUTPUTS.rstfile0 = '''checkpoint'''; +OUTPUTS.job2load = JOB2LOAD; %% Create directories if ~exist(SIMDIR, 'dir') mkdir(SIMDIR) end if ~exist(BASIC.RESDIR, 'dir') mkdir(BASIC.RESDIR) end %% Compile and WRITE input file INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC); nproc = 1; MAKE = 'cd ..; make; cd wk'; system(MAKE); %% disp(['Set up ',SIMID]); disp([resolution,gridname,degngrad]); if RESTART disp(['- restarting from JOBNUM = ',num2str(JOB2LOAD)]); else disp(['- starting from T = 0']); end diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m index 1e2f98d..b5045f9 100644 --- a/matlab/write_fort90.m +++ b/matlab/write_fort90.m @@ -1,76 +1,81 @@ function [INPUT] = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC) % Write the input script "fort.90" with desired parameters INPUT = 'fort.90'; fid = fopen(INPUT,'wt'); fprintf(fid,'&BASIC\n'); fprintf(fid,[' nrun = ', num2str(BASIC.nrun),'\n']); fprintf(fid,[' dt = ', num2str(BASIC.dt),'\n']); fprintf(fid,[' tmax = ', num2str(BASIC.tmax),'\n']); fprintf(fid,[' RESTART = ', BASIC.RESTART,'\n']); fprintf(fid,[' maxruntime = ', num2str(BASIC.maxruntime),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&GRID\n'); fprintf(fid,[' pmaxe =', num2str(GRID.pmaxe),'\n']); fprintf(fid,[' jmaxe = ', num2str(GRID.jmaxe),'\n']); fprintf(fid,[' pmaxi = ', num2str(GRID.pmaxi),'\n']); fprintf(fid,[' jmaxi = ', num2str(GRID.jmaxi),'\n']); fprintf(fid,[' Nr = ', num2str(GRID.Nr),'\n']); fprintf(fid,[' Lr = ', num2str(GRID.Lr),'\n']); fprintf(fid,[' Nz = ', num2str(GRID.Nz),'\n']); fprintf(fid,[' Lz = ', num2str(GRID.Lz),'\n']); fprintf(fid,[' kpar = ', num2str(GRID.kpar),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&OUTPUT_PAR\n'); fprintf(fid,[' nsave_0d = ', num2str(OUTPUTS.nsave_0d),'\n']); fprintf(fid,[' nsave_1d = ', num2str(OUTPUTS.nsave_1d),'\n']); fprintf(fid,[' nsave_2d = ', num2str(OUTPUTS.nsave_2d),'\n']); fprintf(fid,[' nsave_5d = ', num2str(OUTPUTS.nsave_5d),'\n']); fprintf(fid,[' nsave_cp = ', num2str(OUTPUTS.nsave_cp),'\n']); fprintf(fid,[' write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']); +fprintf(fid,[' write_gamma = ', OUTPUTS.write_gamma,'\n']); +fprintf(fid,[' write_phi = ', OUTPUTS.write_phi,'\n']); +fprintf(fid,[' write_Na00 = ', OUTPUTS.write_Na00,'\n']); +fprintf(fid,[' write_Napj = ', OUTPUTS.write_Napj,'\n']); +fprintf(fid,[' write_Sapj = ', OUTPUTS.write_Sapj,'\n']); fprintf(fid,[' resfile0 = ', OUTPUTS.resfile0,'\n']); fprintf(fid,[' rstfile0 = ', OUTPUTS.rstfile0,'\n']); fprintf(fid,[' job2load = ', num2str(OUTPUTS.job2load),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&MODEL_PAR\n'); fprintf(fid,' ! Collisionality\n'); fprintf(fid,[' CO = ', num2str(MODEL.CO),'\n']); fprintf(fid,[' CLOS = ', num2str(MODEL.CLOS),'\n']); fprintf(fid,[' NON_LIN = ', MODEL.NON_LIN,'\n']); fprintf(fid,[' mu = ', num2str(MODEL.mu),'\n']); fprintf(fid,[' mu_p = ', num2str(MODEL.mu_p),'\n']); fprintf(fid,[' mu_j = ', num2str(MODEL.mu_j),'\n']); fprintf(fid,[' nu = ', num2str(MODEL.nu),'\n']); fprintf(fid,[' tau_e = ', num2str(MODEL.tau_e),'\n']); fprintf(fid,[' tau_i = ', num2str(MODEL.tau_i),'\n']); fprintf(fid,[' sigma_e = ', num2str(MODEL.sigma_e),'\n']); fprintf(fid,[' sigma_i = ', num2str(MODEL.sigma_i),'\n']); fprintf(fid,[' q_e = ', num2str(MODEL.q_e),'\n']); fprintf(fid,[' q_i = ', num2str(MODEL.q_i),'\n']); fprintf(fid,[' eta_n = ', num2str(MODEL.eta_n),'\n']); fprintf(fid,[' eta_T = ', num2str(MODEL.eta_T),'\n']); fprintf(fid,[' eta_B = ', num2str(MODEL.eta_B),'\n']); fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']); fprintf(fid,'/\n'); fprintf(fid,'&INITIAL_CON\n'); fprintf(fid,[' INIT_NOISY_PHI =', INITIAL.init_noisy_phi,'\n']); % fprintf(fid,[' only_Na00 =', '.false.','\n']); fprintf(fid,[' initback_moments =', num2str(INITIAL.initback_moments),'\n']); fprintf(fid,[' initnoise_moments =', num2str(INITIAL.initnoise_moments),'\n']); fprintf(fid,[' iseed =', num2str(INITIAL.iseed),'\n']); fprintf(fid,[' selfmat_file =', INITIAL.selfmat_file,'\n']); fprintf(fid,[' eimat_file =', INITIAL.eimat_file,'\n']); fprintf(fid,[' iemat_file =', INITIAL.iemat_file,'\n']); fprintf(fid,'/\n'); fprintf(fid,'&TIME_INTEGRATION_PAR\n'); fprintf(fid,[' numerical_scheme=', TIME_INTEGRATION.numerical_scheme,'\n']); fprintf(fid,'/'); fclose(fid); system(['cp fort.90 ',BASIC.RESDIR,'/.']); end diff --git a/src/basic_mod.F90 b/src/basic_mod.F90 index 15ca635..ecfb416 100644 --- a/src/basic_mod.F90 +++ b/src/basic_mod.F90 @@ -1,291 +1,294 @@ MODULE basic ! Basic module for time dependent problems use, intrinsic :: iso_c_binding use prec_const IMPLICIT none ! INCLUDE 'fftw3-mpi.f03' INTEGER :: nrun = 1 ! Number of time steps to run real(dp) :: tmax = 100000.0 ! Maximum simulation time real(dp) :: dt = 1.0 ! Time step real(dp) :: time = 0 ! Current simulation time (Init from restart file) INTEGER :: comm0 ! Default communicator with a topology INTEGER :: commp, commr ! Communicators for 1-dim cartesian subgrids of comm0 + INTEGER :: commr_p0 ! Communicators along kr for only rank 0 on p INTEGER :: jobnum = 0 ! Job number INTEGER :: step = 0 ! Calculation step of this run INTEGER :: cstep = 0 ! Current step number (Init from restart file) LOGICAL :: RESTART = .FALSE. ! Signal end of run LOGICAL :: nlend = .FALSE. ! Signal end of run LOGICAL :: crashed = .FALSE. ! Signal end of crashed run INTEGER :: ierr ! flag for MPI error INTEGER :: my_id ! Rank in COMM_WORLD INTEGER :: num_procs ! number of MPI processes - INTEGER :: num_procs_p, num_procs_kr ! Number of processes in p and r + INTEGER :: num_procs_p ! Number of processes in p + INTEGER :: num_procs_kr ! Number of processes in r INTEGER :: rank_0, rank_p, rank_r! Ranks in comm0, commp, commr - INTEGER :: nbr_L, nbr_R ! Left and right neighbours (along p) - INTEGER :: nbr_T, nbr_B ! Top and bottom neighbours (along kr) + INTEGER :: nbr_L, nbr_R ! Left and right neighbours (along p) + INTEGER :: nbr_T, nbr_B ! Top and bottom neighbours (along kr) + INTEGER :: iframe0d ! counting the number of times 0d datasets are outputed (for diagnose) INTEGER :: iframe1d ! counting the number of times 1d datasets are outputed (for diagnose) INTEGER :: iframe2d ! counting the number of times 2d datasets are outputed (for diagnose) INTEGER :: iframe3d ! counting the number of times 3d datasets are outputed (for diagnose) INTEGER :: iframe5d ! counting the number of times 5d datasets are outputed (for diagnose) ! List of logical file units INTEGER :: lu_in = 90 ! File duplicated from STDIN INTEGER :: lu_job = 91 ! myjob file ! To measure computation time real :: start, finish real(dp) :: t0_rhs, t0_adv_field, t0_poisson, t0_Sapj, t0_diag, t0_checkfield, t0_step, t0_comm real(dp) :: t1_rhs, t1_adv_field, t1_poisson, t1_Sapj, t1_diag, t1_checkfield, t1_step, t1_comm real(dp) :: tc_rhs, tc_adv_field, tc_poisson, tc_Sapj, tc_diag, tc_checkfield, tc_step, tc_comm real(dp):: maxruntime = 1e9 ! Maximum simulation CPU time INTERFACE allocate_array MODULE PROCEDURE allocate_array_dp1,allocate_array_dp2,allocate_array_dp3,allocate_array_dp4 MODULE PROCEDURE allocate_array_dc1,allocate_array_dc2,allocate_array_dc3,allocate_array_dc4, allocate_array_dc5 MODULE PROCEDURE allocate_array_i1,allocate_array_i2,allocate_array_i3,allocate_array_i4 MODULE PROCEDURE allocate_array_l1,allocate_array_l2,allocate_array_l3,allocate_array_l4 END INTERFACE allocate_array CONTAINS !================================================================================ SUBROUTINE basic_data ! Read basic data for input file use prec_const IMPLICIT NONE NAMELIST /BASIC/ nrun, dt, tmax, RESTART, maxruntime READ(lu_in,basic) !Init cumulative timers tc_rhs = 0.;tc_adv_field = 0.; tc_poisson = 0. tc_Sapj = 0.; tc_diag = 0.; tc_checkfield = 0. END SUBROUTINE basic_data !================================================================================ SUBROUTINE daytim(str) ! Print date and time use prec_const IMPLICIT NONE CHARACTER(len=*) , INTENT(in) :: str CHARACTER(len=16) :: d, t, dat, time !________________________________________________________________________________ ! CALL DATE_AND_TIME(d,t) dat=d(7:8) // '/' // d(5:6) // '/' // d(1:4) time=t(1:2) // ':' // t(3:4) // ':' // t(5:10) WRITE(*,'(a,1x,a,1x,a)') str, dat(1:10), time(1:12) ! END SUBROUTINE daytim !================================================================================ SUBROUTINE display_h_min_s(time) real :: time integer :: days, hours, mins, secs days = FLOOR(time/24./3600.); hours= FLOOR(time/3600.); mins = FLOOR(time/60.); secs = FLOOR(time); IF ( days .GT. 0 ) THEN !display day h min s hours = (time/3600./24. - days) * 24 mins = (time/3600. - days*24. - hours) * 60 secs = (time/60. - days*24.*60 - hours*60 - mins) * 60 WRITE(*,*) 'CPU Time = ', days, '[day]', hours, '[h]', mins, '[min]', secs, '[s]' WRITE(*,*) '(',time,'[s])' ELSEIF ( hours .GT. 0 ) THEN !display h min s mins = (time/3600. - hours) * 60 secs = (time/60. - hours*60 - mins) * 60 WRITE(*,*) 'CPU Time = ', hours, '[h]', mins, '[min]', secs, '[s]' WRITE(*,*) '(',time,'[s])' ELSEIF ( mins .GT. 0 ) THEN !display min s secs = (time/60. - mins) * 60 WRITE(*,*) 'CPU Time = ', mins, '[min]', secs, '[s]' WRITE(*,*) '(',time,'[s])' ELSE ! display s WRITE(*,*) 'CPU Time = ', FLOOR(time), '[s]' ENDIF END SUBROUTINE display_h_min_s !================================================================================ ! To allocate arrays of doubles, integers, etc. at run time SUBROUTINE allocate_array_dp1(a,is1,ie1) IMPLICIT NONE real(dp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1 ALLOCATE(a(is1:ie1)) a=0.0_dp END SUBROUTINE allocate_array_dp1 SUBROUTINE allocate_array_dp2(a,is1,ie1,is2,ie2) IMPLICIT NONE real(dp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2 ALLOCATE(a(is1:ie1,is2:ie2)) a=0.0_dp END SUBROUTINE allocate_array_dp2 SUBROUTINE allocate_array_dp3(a,is1,ie1,is2,ie2,is3,ie3) IMPLICIT NONE real(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3)) a=0.0_dp END SUBROUTINE allocate_array_dp3 SUBROUTINE allocate_array_dp4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4) IMPLICIT NONE real(dp), DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4)) a=0.0_dp END SUBROUTINE allocate_array_dp4 SUBROUTINE allocate_array_dp5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5) IMPLICIT NONE real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5)) a=0.0_dp END SUBROUTINE allocate_array_dp5 !======================================== SUBROUTINE allocate_array_dc1(a,is1,ie1) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1 ALLOCATE(a(is1:ie1)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc1 SUBROUTINE allocate_array_dc2(a,is1,ie1,is2,ie2) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2 ALLOCATE(a(is1:ie1,is2:ie2)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc2 SUBROUTINE allocate_array_dc3(a,is1,ie1,is2,ie2,is3,ie3) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc3 SUBROUTINE allocate_array_dc4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc4 SUBROUTINE allocate_array_dc5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc5 !======================================== SUBROUTINE allocate_array_i1(a,is1,ie1) IMPLICIT NONE INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1 ALLOCATE(a(is1:ie1)) a=0 END SUBROUTINE allocate_array_i1 SUBROUTINE allocate_array_i2(a,is1,ie1,is2,ie2) IMPLICIT NONE INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2 ALLOCATE(a(is1:ie1,is2:ie2)) a=0 END SUBROUTINE allocate_array_i2 SUBROUTINE allocate_array_i3(a,is1,ie1,is2,ie2,is3,ie3) IMPLICIT NONE INTEGER, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3)) a=0 END SUBROUTINE allocate_array_i3 SUBROUTINE allocate_array_i4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4) IMPLICIT NONE INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4)) a=0 END SUBROUTINE allocate_array_i4 SUBROUTINE allocate_array_i5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5) IMPLICIT NONE real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5)) a=0 END SUBROUTINE allocate_array_i5 !======================================== SUBROUTINE allocate_array_l1(a,is1,ie1) IMPLICIT NONE LOGICAL, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1 ALLOCATE(a(is1:ie1)) a=.false. END SUBROUTINE allocate_array_l1 SUBROUTINE allocate_array_l2(a,is1,ie1,is2,ie2) IMPLICIT NONE LOGICAL, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2 ALLOCATE(a(is1:ie1,is2:ie2)) a=.false. END SUBROUTINE allocate_array_l2 SUBROUTINE allocate_array_l3(a,is1,ie1,is2,ie2,is3,ie3) IMPLICIT NONE LOGICAL, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3)) a=.false. END SUBROUTINE allocate_array_l3 SUBROUTINE allocate_array_l4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4) IMPLICIT NONE LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4)) a=.false. END SUBROUTINE allocate_array_l4 SUBROUTINE allocate_array_l5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5) IMPLICIT NONE real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5)) a=.false. END SUBROUTINE allocate_array_l5 END MODULE basic diff --git a/src/diagnose.F90 b/src/diagnose.F90 index f4d9a28..40094b6 100644 --- a/src/diagnose.F90 +++ b/src/diagnose.F90 @@ -1,468 +1,519 @@ SUBROUTINE diagnose(kstep) ! Diagnostics, writing simulation state to disk USE basic USE grid USE diagnostics_par USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf, putarrnd USE model USE initial_par USE fields USE time_integration USE prec_const IMPLICIT NONE INCLUDE 'srcinfo.h' INTEGER, INTENT(in) :: kstep INTEGER, parameter :: BUFSIZE = 2 - INTEGER :: rank, dims(1) = (/0/) + INTEGER :: rank = 0 + INTEGER :: dims(1) = (/0/) INTEGER :: cp_counter = 0 CHARACTER(len=256) :: str, fname,test_ ! putarr(...,pardim=1) does not work for 2D domain decomposition ! so we need to gather non 5D data on one proc to output it INTEGER :: parray_e_full(1:pmaxe+1), parray_i_full(1:pmaxi+1) INTEGER :: jarray_e_full(1:jmaxe+1), jarray_i_full(1:jmaxi+1) REAL(dp) :: krarray_full(1:nkr), kzarray_full(1:nkz) CALL cpu_time(t0_diag) ! Measuring time !_____________________________________________________________________________ ! 1. Initial diagnostics IF ((kstep .EQ. 0)) THEN ! Writing output filename WRITE(resfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',jobnum,'.h5' ! 1.1 Initial run ! Main output file creation IF (write_doubleprecision) THEN CALL creatf(resfile, fidres, real_prec='d', mpicomm=comm0) ELSE CALL creatf(resfile, fidres, mpicomm=comm0) END IF IF (my_id .EQ. 0) WRITE(*,'(3x,a,a)') TRIM(resfile), ' created' ! Data group CALL creatg(fidres, "/data", "data") - CALL creatg(fidres, "/data/var2d", "2d profiles") - CALL creatg(fidres, "/data/var5d", "5d profiles") - - ! Initialize counter of number of saves for each category - IF (cstep==0) THEN - iframe2d=0 - ENDIF - CALL attach(fidres,"/data/var2d/" , "frames", iframe2d) - IF (cstep==0) THEN - iframe5d=0 - END IF - CALL attach(fidres,"/data/var5d/" , "frames", iframe5d) ! File group CALL creatg(fidres, "/files", "files") CALL attach(fidres, "/files", "jobnum", jobnum) ! Profiler time measurement CALL creatg(fidres, "/profiler", "performance analysis") CALL creatd(fidres, 0, dims, "/profiler/Tc_rhs", "cumulative rhs computation time") CALL creatd(fidres, 0, dims, "/profiler/Tc_adv_field", "cumulative adv. fields computation time") CALL creatd(fidres, 0, dims, "/profiler/Tc_comm", "cumulative communication time") CALL creatd(fidres, 0, dims, "/profiler/Tc_poisson", "cumulative poisson computation time") CALL creatd(fidres, 0, dims, "/profiler/Tc_Sapj", "cumulative Sapj computation time") CALL creatd(fidres, 0, dims, "/profiler/Tc_diag", "cumulative sym computation time") CALL creatd(fidres, 0, dims, "/profiler/Tc_checkfield", "cumulative checkfield computation time") CALL creatd(fidres, 0, dims, "/profiler/Tc_step", "cumulative total step computation time") CALL creatd(fidres, 0, dims, "/profiler/time", "current simulation time") ! Build the full grids on process 0 to diagnose it without comm IF (my_id .EQ. 0) THEN ! P DO ip = 1,pmaxe+1; parray_e_full(ip) = (ip-1); END DO DO ip = 1,pmaxi+1; parray_i_full(ip) = (ip-1); END DO ! J DO ij = 1,jmaxe+1; jarray_e_full(ij) = (ij-1); END DO DO ij = 1,jmaxi+1; jarray_i_full(ij) = (ij-1); END DO ! Kr DO ikr = 1,Nkr krarray_full(ikr) = REAL(ikr-1,dp) * deltakr END DO ! Kz IF (Nkz .GT. 1) THEN DO ikz = 1,Nkz kzarray_full(ikz) = deltakz*(MODULO(ikz-1,Nkz/2)-Nkz/2*FLOOR(2.*real(ikz-1)/real(Nkz))) if (ikz .EQ. Nz/2+1) kzarray(ikz) = -kzarray(ikz) END DO ELSE kzarray_full(1) = 0 endif ENDIF - ! var2d group (electro. pot., Ni00 moment) - rank = 0 - CALL creatd(fidres, rank, dims, "/data/var2d/time", "Time t*c_s/R") - CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number") + ! Grid info CALL creatg(fidres, "/data/grid", "Grid data") CALL putarr(fidres, "/data/grid/coordkr", krarray_full(1:nkr),"kr*rho_s0", ionode=0) CALL putarr(fidres, "/data/grid/coordkz", kzarray_full(1:nkz),"kz*rho_s0", ionode=0) - CALL putarr(fidres, "/data/grid/coordp" , parray_e_full(1:pmaxe+1), "p_e", ionode=0) - CALL putarr(fidres, "/data/grid/coordj" , jarray_e_full(1:jmaxe+1), "j_e", ionode=0) + CALL putarr(fidres, "/data/grid/coordp_e" , parray_e_full(1:pmaxe+1), "p_e", ionode=0) + CALL putarr(fidres, "/data/grid/coordj_e" , jarray_e_full(1:jmaxe+1), "j_e", ionode=0) + CALL putarr(fidres, "/data/grid/coordp_i" , parray_i_full(1:pmaxe+1), "p_i", ionode=0) + CALL putarr(fidres, "/data/grid/coordj_i" , jarray_i_full(1:jmaxe+1), "j_i", ionode=0) + + ! var0d group (gyro transport) + IF (nsave_0d .GT. 0) THEN + CALL creatg(fidres, "/data/var0d", "0d profiles") + CALL creatd(fidres, rank, dims, "/data/var0d/time", "Time t*c_s/R") + CALL creatd(fidres, rank, dims, "/data/var0d/cstep", "iteration number") + + IF (write_gamma) THEN + CALL creatd(fidres, rank, dims, "/data/var0d/gflux_ri", "Radial gyro ion transport") + CALL creatd(fidres, rank, dims, "/data/var0d/pflux_ri", "Radial part ion transport") + ENDIF + IF (cstep==0) THEN + iframe0d=0 + ENDIF + CALL attach(fidres,"/data/var0d/" , "frames", iframe0d) + END IF + + ! var2d group (electro. pot., Ni00 moment) IF (nsave_2d .GT. 0) THEN + CALL creatg(fidres, "/data/var2d", "2d profiles") + CALL creatd(fidres, rank, dims, "/data/var2d/time", "Time t*c_s/R") + CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number") + + IF (write_phi) CALL creatg(fidres, "/data/var2d/phi", "phi") + + IF (write_Na00) THEN CALL creatg(fidres, "/data/var2d/Ne00", "Ne00") CALL creatg(fidres, "/data/var2d/Ni00", "Ni00") - CALL creatg(fidres, "/data/var2d/phi", "phi") + ENDIF + + IF (cstep==0) THEN + iframe2d=0 + ENDIF + CALL attach(fidres,"/data/var2d/" , "frames", iframe2d) END IF ! var5d group (moments) - rank = 0 - CALL creatd(fidres, rank, dims, "/data/var5d/time", "Time t*c_s/R") - CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number") IF (nsave_5d .GT. 0) THEN - CALL creatg(fidres, "/data/var5d/moments_e", "moments_e") - CALL creatg(fidres, "/data/var5d/moments_i", "moments_i") - ! CALL creatg(fidres, "/data/var5d/Sepj", "Sepj") - ! CALL creatg(fidres, "/data/var5d/Sipj", "Sipj") + CALL creatg(fidres, "/data/var5d", "5d profiles") + CALL creatd(fidres, rank, dims, "/data/var5d/time", "Time t*c_s/R") + CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number") + + IF (write_Napj) THEN + CALL creatg(fidres, "/data/var5d/moments_e", "moments_e") + CALL creatg(fidres, "/data/var5d/moments_i", "moments_i") + ENDIF + + IF (write_Sapj) THEN + CALL creatg(fidres, "/data/var5d/moments_e", "Sipj") + CALL creatg(fidres, "/data/var5d/moments_i", "Sepj") + ENDIF + + IF (cstep==0) THEN + iframe5d=0 + END IF + CALL attach(fidres,"/data/var5d/" , "frames", iframe5d) END IF ! Add input namelist variables as attributes of /data/input, defined in srcinfo.h IF (my_id .EQ. 0) WRITE(*,*) 'VERSION=', VERSION IF (my_id .EQ. 0) WRITE(*,*) 'BRANCH=', BRANCH IF (my_id .EQ. 0) WRITE(*,*) 'AUTHOR=', AUTHOR IF (my_id .EQ. 0) WRITE(*,*) 'HOST=', HOST WRITE(str,'(a,i2.2)') "/data/input" - - rank=0 CALL creatd(fidres, rank,dims,TRIM(str),'Input parameters') CALL attach(fidres, TRIM(str), "version", VERSION) !defined in srcinfo.h CALL attach(fidres, TRIM(str), "branch", BRANCH) !defined in srcinfo.h CALL attach(fidres, TRIM(str), "author", AUTHOR) !defined in srcinfo.h CALL attach(fidres, TRIM(str), "execdate", EXECDATE) !defined in srcinfo.h CALL attach(fidres, TRIM(str), "host", HOST) !defined in srcinfo.h CALL attach(fidres, TRIM(str), "start_time", time) CALL attach(fidres, TRIM(str), "start_cstep", cstep-1) + CALL attach(fidres, TRIM(str), "start_iframe0d", iframe0d) CALL attach(fidres, TRIM(str), "start_iframe2d", iframe2d) CALL attach(fidres, TRIM(str), "start_iframe5d", iframe5d) CALL attach(fidres, TRIM(str), "dt", dt) CALL attach(fidres, TRIM(str), "tmax", tmax) CALL attach(fidres, TRIM(str), "nrun", nrun) CALL attach(fidres, TRIM(str), "cpu_time", -1) CALL attach(fidres, TRIM(str), "Nproc", num_procs) CALL attach(fidres, TRIM(str), "Np_p" , num_procs_p) CALL attach(fidres, TRIM(str), "Np_kr",num_procs_kr) + CALL attach(fidres, TRIM(str), "write_gamma",write_gamma) + CALL attach(fidres, TRIM(str), "write_phi",write_phi) + CALL attach(fidres, TRIM(str), "write_Na00",write_Na00) + CALL attach(fidres, TRIM(str), "write_Napj",write_Napj) + CALL attach(fidres, TRIM(str), "write_Sapj",write_Sapj) CALL grid_outputinputs(fidres, str) CALL output_par_outputinputs(fidres, str) CALL model_outputinputs(fidres, str) CALL initial_outputinputs(fidres, str) CALL time_integration_outputinputs(fidres, str) ! Save STDIN (input file) of this run IF(jobnum .LE. 99) THEN WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum ELSE WRITE(str,'(a,i3.2)') "/files/STDIN.",jobnum END IF INQUIRE(unit=lu_in, name=fname) CLOSE(lu_in) CALL putfile(fidres, TRIM(str), TRIM(fname),ionode=0) ELSEIF((kstep .EQ. 0)) THEN IF(jobnum .LE. 99) THEN WRITE(resfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',jobnum,'.h5' ELSE WRITE(resfile,'(a,a1,i3.2,a3)') TRIM(resfile0),'_',jobnum,'.h5' END IF CALL openf(resfile,fidres, 'D'); ENDIF !_____________________________________________________________________________ ! 2. Periodic diagnostics ! IF (kstep .GE. 0) THEN ! Terminal info IF (MOD(cstep, INT(1.0/dt)) == 0 .AND. (my_id .EQ. 0)) THEN WRITE(*,"(F5.0,A,F5.0)") time,"/",tmax ENDIF ! 2.1 0d history arrays IF (nsave_0d .GT. 0) THEN IF ( MOD(cstep, nsave_0d) == 0 ) THEN CALL diagnose_0d END IF END IF ! 2.2 1d profiles ! empty in our case ! 2.3 2d profiles IF (nsave_2d .GT. 0) THEN IF (MOD(cstep, nsave_2d) == 0) THEN CALL diagnose_2d END IF END IF ! 2.4 3d profiles IF (nsave_5d .GT. 0 .AND. cstep .GT. 0) THEN IF (MOD(cstep, nsave_5d) == 0) THEN CALL diagnose_5d END IF END IF !_____________________________________________________________________________ ! 3. Final diagnostics ELSEIF (kstep .EQ. -1) THEN CALL cpu_time(finish) CALL attach(fidres, "/data/input","cpu_time",finish-start) ! Display computational time cost IF (my_id .EQ. 0) CALL display_h_min_s(finish-start) ! Close all diagnostic files CALL closef(fidres) END IF CALL cpu_time(t1_diag); tc_diag = tc_diag + (t1_diag - t0_diag) END SUBROUTINE diagnose SUBROUTINE diagnose_0d USE basic - USE futils, ONLY: append + USE futils, ONLY: append, attach, getatt USE diagnostics_par USE prec_const + USE processing IMPLICIT NONE - + ! Time measurement data CALL append(fidres, "/profiler/Tc_rhs", tc_rhs,ionode=0) CALL append(fidres, "/profiler/Tc_adv_field", tc_adv_field,ionode=0) CALL append(fidres, "/profiler/Tc_poisson", tc_poisson,ionode=0) CALL append(fidres, "/profiler/Tc_Sapj", tc_Sapj,ionode=0) CALL append(fidres, "/profiler/Tc_diag", tc_diag,ionode=0) CALL append(fidres, "/profiler/Tc_checkfield",tc_checkfield,ionode=0) CALL append(fidres, "/profiler/Tc_comm", tc_comm,ionode=0) CALL append(fidres, "/profiler/Tc_step", tc_step,ionode=0) CALL append(fidres, "/profiler/time", time,ionode=0) + ! Processing data + CALL append(fidres, "/data/var0d/time", time,ionode=0) + CALL append(fidres, "/data/var0d/cstep", real(cstep,dp),ionode=0) + CALL getatt(fidres, "/data/var0d/", "frames",iframe2d) + iframe0d=iframe0d+1 + CALL attach(fidres,"/data/var0d/" , "frames", iframe0d) + ! Ion transport data + IF (write_gamma) THEN + CALL compute_radial_ion_transport + CALL append(fidres, "/data/var0d/gflux_ri",gflux_ri,ionode=0) + CALL append(fidres, "/data/var0d/pflux_ri",pflux_ri,ionode=0) + ENDIF END SUBROUTINE diagnose_0d SUBROUTINE diagnose_2d USE basic USE futils, ONLY: append, getatt, attach, putarrnd USE fields USE array, ONLY: Ne00, Ni00 USE grid, ONLY: ikrs,ikre, ikzs,ikze, nkr, nkz, local_nkr, ikr, ikz, ips_e, ips_i USE time_integration USE diagnostics_par USE prec_const IMPLICIT NONE COMPLEX(dp) :: buffer(ikrs:ikre,ikzs:ikze) INTEGER :: i_, root, world_rank, world_size CALL append(fidres, "/data/var2d/time", time,ionode=0) CALL append(fidres, "/data/var2d/cstep", real(cstep,dp),ionode=0) CALL getatt(fidres, "/data/var2d/", "frames",iframe2d) iframe2d=iframe2d+1 CALL attach(fidres,"/data/var2d/" , "frames", iframe2d) - CALL write_field2d(phi (:,:), 'phi') + IF (write_phi) CALL write_field2d(phi (:,:), 'phi') - IF ( (ips_e .EQ. 1) .AND. (ips_i .EQ. 1) ) THEN - Ne00(ikrs:ikre,ikzs:ikze) = moments_e(ips_e,1,ikrs:ikre,ikzs:ikze,updatetlevel) - Ni00(ikrs:ikre,ikzs:ikze) = moments_i(ips_e,1,ikrs:ikre,ikzs:ikze,updatetlevel) - ENDIF - - root = 0 - - !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!! - CALL MPI_COMM_RANK(commp,world_rank,ierr) - CALL MPI_COMM_SIZE(commp,world_size,ierr) + IF (write_Na00) THEN + IF ( (ips_e .EQ. 1) .AND. (ips_i .EQ. 1) ) THEN + Ne00(ikrs:ikre,ikzs:ikze) = moments_e(ips_e,1,ikrs:ikre,ikzs:ikze,updatetlevel) + Ni00(ikrs:ikre,ikzs:ikze) = moments_i(ips_e,1,ikrs:ikre,ikzs:ikze,updatetlevel) + ENDIF - IF (world_size .GT. 1) THEN - !! Broadcast phi to the other processes on the same k range (communicator along p) - IF (world_rank .EQ. root) THEN - ! Fill the buffer - DO ikr = ikrs,ikre - DO ikz = ikzs,ikze - buffer(ikr,ikz) = Ne00(ikr,ikz) + root = 0 + !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!! + CALL MPI_COMM_RANK(commp,world_rank,ierr) + CALL MPI_COMM_SIZE(commp,world_size,ierr) + + IF (world_size .GT. 1) THEN + !! Broadcast phi to the other processes on the same k range (communicator along p) + IF (world_rank .EQ. root) THEN + ! Fill the buffer + DO ikr = ikrs,ikre + DO ikz = ikzs,ikze + buffer(ikr,ikz) = Ne00(ikr,ikz) + ENDDO ENDDO - ENDDO - ! Send it to all the other processes - DO i_ = 0,num_procs_p-1 - IF (i_ .NE. world_rank) & - CALL MPI_SEND(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, i_, 0, commp, ierr) - ENDDO - ELSE - ! Recieve buffer from root - CALL MPI_RECV(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, root, 0, commp, MPI_STATUS_IGNORE, ierr) - ! Write it in phi - DO ikr = ikrs,ikre - DO ikz = ikzs,ikze - Ne00(ikr,ikz) = buffer(ikr,ikz) + ! Send it to all the other processes + DO i_ = 0,num_procs_p-1 + IF (i_ .NE. world_rank) & + CALL MPI_SEND(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, i_, 0, commp, ierr) + ENDDO + ELSE + ! Recieve buffer from root + CALL MPI_RECV(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, root, 0, commp, MPI_STATUS_IGNORE, ierr) + ! Write it in phi + DO ikr = ikrs,ikre + DO ikz = ikzs,ikze + Ne00(ikr,ikz) = buffer(ikr,ikz) + ENDDO ENDDO - ENDDO + ENDIF ENDIF - ENDIF - CALL write_field2d(Ne00(ikrs:ikre,ikzs:ikze), 'Ne00') + CALL write_field2d(Ne00(ikrs:ikre,ikzs:ikze), 'Ne00') - !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!! - CALL MPI_COMM_RANK(commp,world_rank,ierr) - CALL MPI_COMM_SIZE(commp,world_size,ierr) - - IF (world_size .GT. 1) THEN - !! Broadcast phi to the other processes on the same k range (communicator along p) - IF (world_rank .EQ. root) THEN - ! Fill the buffer - DO ikr = ikrs,ikre - DO ikz = ikzs,ikze - buffer(ikr,ikz) = Ni00(ikr,ikz) + !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!! + CALL MPI_COMM_RANK(commp,world_rank,ierr) + CALL MPI_COMM_SIZE(commp,world_size,ierr) + + IF (world_size .GT. 1) THEN + !! Broadcast phi to the other processes on the same k range (communicator along p) + IF (world_rank .EQ. root) THEN + ! Fill the buffer + DO ikr = ikrs,ikre + DO ikz = ikzs,ikze + buffer(ikr,ikz) = Ni00(ikr,ikz) + ENDDO ENDDO - ENDDO - ! Send it to all the other processes - DO i_ = 0,num_procs_p-1 - IF (i_ .NE. world_rank) & - CALL MPI_SEND(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, i_, 0, commp, ierr) - ENDDO - ELSE - ! Recieve buffer from root - CALL MPI_RECV(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, root, 0, commp, MPI_STATUS_IGNORE, ierr) - ! Write it in phi - DO ikr = ikrs,ikre - DO ikz = ikzs,ikze - Ni00(ikr,ikz) = buffer(ikr,ikz) + ! Send it to all the other processes + DO i_ = 0,num_procs_p-1 + IF (i_ .NE. world_rank) & + CALL MPI_SEND(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, i_, 0, commp, ierr) + ENDDO + ELSE + ! Recieve buffer from root + CALL MPI_RECV(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, root, 0, commp, MPI_STATUS_IGNORE, ierr) + ! Write it in phi + DO ikr = ikrs,ikre + DO ikz = ikzs,ikze + Ni00(ikr,ikz) = buffer(ikr,ikz) + ENDDO ENDDO - ENDDO + ENDIF ENDIF - ENDIF - CALL write_field2d(Ni00(ikrs:ikre,ikzs:ikze), 'Ni00') + CALL write_field2d(Ni00(ikrs:ikre,ikzs:ikze), 'Ni00') + ENDIF CONTAINS SUBROUTINE write_field2d(field, text) USE futils, ONLY: attach, putarr USE grid, ONLY: ikrs,ikre, ikzs,ikze, nkr, nkz, local_nkr USE prec_const USE basic, ONLY : commr, num_procs_p, rank_p IMPLICIT NONE COMPLEX(dp), DIMENSION(ikrs:ikre, ikzs:ikze), INTENT(IN) :: field CHARACTER(*), INTENT(IN) :: text COMPLEX(dp) :: buffer_dist(ikrs:ikre,ikzs:ikze) COMPLEX(dp) :: buffer_full(1:nkr,1:nkz) INTEGER :: scount, rcount CHARACTER(LEN=50) :: dset_name scount = (ikre-ikrs+1) * (ikze-ikzs+1) rcount = scount WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var2d", TRIM(text), iframe2d IF (num_procs .EQ. 1) THEN ! no data distribution CALL putarr(fidres, dset_name, field(ikrs:ikre, ikzs:ikze), ionode=0) ELSE CALL putarrnd(fidres, dset_name, field(ikrs:ikre, ikzs:ikze), (/1, 1/)) ENDIF CALL attach(fidres, dset_name, "time", time) END SUBROUTINE write_field2d END SUBROUTINE diagnose_2d SUBROUTINE diagnose_5d USE basic USE futils, ONLY: append, getatt, attach, putarrnd USE fields USE array!, ONLY: Sepj, Sipj USE grid, ONLY: ips_e,ipe_e, ips_i, ipe_i, & ijs_e,ije_e, ijs_i, ije_i USE time_integration USE diagnostics_par USE prec_const IMPLICIT NONE CALL append(fidres, "/data/var5d/time", time,ionode=0) CALL append(fidres, "/data/var5d/cstep", real(cstep,dp),ionode=0) CALL getatt(fidres, "/data/var5d/", "frames",iframe5d) iframe5d=iframe5d+1 CALL attach(fidres,"/data/var5d/" , "frames", iframe5d) - CALL write_field5d_e(moments_e(ips_e:ipe_e,ijs_e:ije_e,:,:,updatetlevel), 'moments_e') - CALL write_field5d_i(moments_i(ips_i:ipe_i,ijs_i:ije_i,:,:,updatetlevel), 'moments_i') + IF (write_Napj) THEN + CALL write_field5d_e(moments_e(ips_e:ipe_e,ijs_e:ije_e,:,:,updatetlevel), 'moments_e') + CALL write_field5d_i(moments_i(ips_i:ipe_i,ijs_i:ije_i,:,:,updatetlevel), 'moments_i') + ENDIF - ! CALL write_field5d_e(Sepj(ips_e:ipe_e,ijs_e:ije_e,:,:), 'Sepj') - ! CALL write_field5d_i(Sipj(ips_i:ipe_i,ijs_i:ije_i,:,:), 'Sipj') + IF (write_Sapj) THEN + CALL write_field5d_e(Sepj(ips_e:ipe_e,ijs_e:ije_e,:,:), 'Sepj') + CALL write_field5d_i(Sipj(ips_i:ipe_i,ijs_i:ije_i,:,:), 'Sipj') + ENDIF CONTAINS SUBROUTINE write_field5d_e(field, text) USE futils, ONLY: attach, putarr, putarrnd USE grid, ONLY: ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze USE prec_const IMPLICIT NONE COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze), INTENT(IN) :: field CHARACTER(*), INTENT(IN) :: text CHARACTER(LEN=50) :: dset_name WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d IF (num_procs .EQ. 1) THEN CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze), ionode=0) ELSE CALL putarrnd(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze), (/1,3/)) ENDIF CALL attach(fidres, dset_name, 'cstep', cstep) CALL attach(fidres, dset_name, 'time', time) CALL attach(fidres, dset_name, 'jobnum', jobnum) CALL attach(fidres, dset_name, 'dt', dt) CALL attach(fidres, dset_name, 'iframe2d', iframe2d) CALL attach(fidres, dset_name, 'iframe5d', iframe5d) END SUBROUTINE write_field5d_e SUBROUTINE write_field5d_i(field, text) USE futils, ONLY: attach, putarr, putarrnd USE grid, ONLY: ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze USE prec_const IMPLICIT NONE COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze), INTENT(IN) :: field CHARACTER(*), INTENT(IN) :: text CHARACTER(LEN=50) :: dset_name WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d IF (num_procs .EQ. 1) THEN CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze), ionode=0) ELSE CALL putarrnd(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze), (/1,3/)) ENDIF CALL attach(fidres, dset_name, 'cstep', cstep) CALL attach(fidres, dset_name, 'time', time) CALL attach(fidres, dset_name, 'jobnum', jobnum) CALL attach(fidres, dset_name, 'dt', dt) CALL attach(fidres, dset_name, 'iframe2d', iframe2d) CALL attach(fidres, dset_name, 'iframe5d', iframe5d) END SUBROUTINE write_field5d_i END SUBROUTINE diagnose_5d diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90 index a0fb7a6..0f74389 100644 --- a/src/diagnostics_par_mod.F90 +++ b/src/diagnostics_par_mod.F90 @@ -1,62 +1,66 @@ MODULE diagnostics_par ! Module for diagnostic parameters USE prec_const IMPLICIT NONE PRIVATE - LOGICAL, PUBLIC, PROTECTED :: write_doubleprecision=.FALSE. + LOGICAL, PUBLIC, PROTECTED :: write_doubleprecision = .FALSE. + LOGICAL, PUBLIC, PROTECTED :: write_gamma + LOGICAL, PUBLIC, PROTECTED :: write_phi, write_Na00 + LOGICAL, PUBLIC, PROTECTED :: write_Napj, write_Sapj INTEGER, PUBLIC, PROTECTED :: nsave_0d , nsave_1d , nsave_2d , nsave_5d, nsave_cp ! HDF5 file CHARACTER(len=256), PUBLIC :: resfile0 = "results" ! Head of main result file name CHARACTER(len=256), PUBLIC :: resfile ! Main result file INTEGER, PUBLIC :: job2load ! jobnum of the checkpoint to load INTEGER, PUBLIC :: fidres ! FID for resfile CHARACTER(len=256), PUBLIC :: rstfile0 = "restart" ! Head of restart file name CHARACTER(len=256), PUBLIC :: rstfile ! Full restart file INTEGER, PUBLIC :: fidrst ! FID for restart file PUBLIC :: output_par_readinputs, output_par_outputinputs CONTAINS SUBROUTINE output_par_readinputs ! Read the input parameters USE basic, ONLY : lu_in USE prec_const IMPLICIT NONE NAMELIST /OUTPUT_PAR/ nsave_0d , nsave_1d , nsave_2d , nsave_5d, nsave_cp - NAMELIST /OUTPUT_PAR/ write_doubleprecision + NAMELIST /OUTPUT_PAR/ write_doubleprecision, write_gamma, write_phi + NAMELIST /OUTPUT_PAR/ write_Na00, write_Napj, write_Sapj NAMELIST /OUTPUT_PAR/ resfile0, rstfile0, job2load READ(lu_in,output_par) END SUBROUTINE output_par_readinputs SUBROUTINE output_par_outputinputs(fidres, str) ! ! Write the input parameters to the results_xx.h5 file ! USE prec_const USE futils, ONLY: attach IMPLICIT NONE INTEGER, INTENT(in) :: fidres CHARACTER(len=256), INTENT(in) :: str CALL attach(fidres, TRIM(str), "write_doubleprecision", write_doubleprecision) CALL attach(fidres, TRIM(str), "nsave_0d", nsave_0d) CALL attach(fidres, TRIM(str), "nsave_1d", nsave_1d) CALL attach(fidres, TRIM(str), "nsave_2d", nsave_2d) CALL attach(fidres, TRIM(str), "nsave_5d", nsave_5d) CALL attach(fidres, TRIM(str), "resfile0", resfile0) END SUBROUTINE output_par_outputinputs END MODULE diagnostics_par diff --git a/src/processing_mod.F90 b/src/processing_mod.F90 new file mode 100644 index 0000000..a99c567 --- /dev/null +++ b/src/processing_mod.F90 @@ -0,0 +1,65 @@ +MODULE processing + ! contains the Hermite-Laguerre collision operators. Solved using COSOlver. + USE basic + USE prec_const + USE grid + implicit none + + REAL(dp), PUBLIC, PROTECTED :: pflux_ri, gflux_ri + + PUBLIC :: compute_radial_ion_transport + +CONTAINS + +SUBROUTINE compute_radial_ion_transport + USE fields, ONLY : moments_i, phi + USE array, ONLY : kernel_i + USE time_integration, ONLY : updatetlevel + IMPLICIT NONE + COMPLEX(dp) :: pflux_local, gflux_local + REAL(dp) :: kz_, buffer(1:2) + INTEGER :: i_, world_rank, world_size, root + + pflux_local = 0._dp + gflux_local = 0._dp + IF(ips_i .EQ. 1) THEN + ! Loop to compute gamma_kr = sum_kz sum_j -i k_z Kernel_j Ni00 * phi + DO ikz = ikzs,ikze + kz_ = kzarray(ikz) + DO ikr = ikrs,ikre + gflux_local = gflux_local - & + imagu * kz_ * moments_i(1,1,ikr,ikz,updatetlevel) * CONJG(phi(ikr,ikz)) + DO ij = ijs_i, ije_i + pflux_local = pflux_local - & + imagu * kz_ * kernel_i(ij,ikr,ikz) * moments_i(1,ij,ikr,ikz,updatetlevel) * CONJG(phi(ikr,ikz)) + ENDDO + ENDDO + ENDDO + + buffer(1) = REAL(gflux_local) + buffer(2) = REAL(pflux_local) + + root = 0 + !Gather manually among the rank_p=0 processes and perform the sum + gflux_ri = 0 + pflux_ri = 0 + IF (num_procs_kr .GT. 1) THEN + !! Everyone sends its local_sum to root = 0 + IF (rank_r .NE. root) THEN + CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, commr, ierr) + ELSE + ! Recieve from all the other processes + DO i_ = 0,num_procs_kr-1 + IF (i_ .NE. rank_r) & + CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, commr, MPI_STATUS_IGNORE, ierr) + gflux_ri = gflux_ri + buffer(1) + pflux_ri = pflux_ri + buffer(2) + ENDDO + ENDIF + ENDIF + ENDIF + + +END SUBROUTINE compute_radial_ion_transport + +END MODULE processing \ No newline at end of file diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m index a9b1226..bd544c9 100644 --- a/wk/analysis_2D.m +++ b/wk/analysis_2D.m @@ -1,408 +1,415 @@ %% Load results outfile =''; if 0 %% Load from Marconi outfile =''; outfile =''; outfile =''; outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.4_2_12_eta_0.6_nu_1e-01/50x25_L_100_P_20_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_3e-01/out.txt'; BASIC.RESDIR = load_marconi(outfile); end if 0 %% Load from Daint outfile =''; BASIC.RESDIR = load_daint(outfile); end %% % JOBNUM = 0; load_results; % JOBNUM = 1; load_results; compile_results load_params %% 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 Nkr = numel(kr); Nkz = numel(kz); [KZ,KR] = meshgrid(kz,kr); Lkr = max(kr)-min(kr); Lkz = max(kz)-min(kz); dkr = Lkr/(Nkr-1); dkz = Lkz/(Nkz-1); KPERP2 = KZ.^2+KR.^2; [~,ikr0] = min(abs(kr)); [~,ikz0] = min(abs(kz)); Lk = max(Lkr,Lkz); dr = 2*pi/Lk; dz = 2*pi/Lk; Nr = max(Nkr,Nkz); Nz = Nr; r = dr*(-Nr/2:(Nr/2-1)); Lr = max(r)-min(r); z = dz*(-Nz/2:(Nz/2-1)); Lz = 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(Nr,Nz,Ns2D); % Gyrocenter density ni00 = zeros(Nr,Nz,Ns2D); np_i = zeros(Nr,Nz,Ns5D); % Ion particle density si00 = zeros(Nr,Nz,Ns5D); phi = zeros(Nr,Nz,Ns2D); drphi = zeros(Nr,Nz,Ns2D); dr2phi = zeros(Nr,Nz,Ns2D); dzphi = zeros(Nr,Nz,Ns2D); for it = 1:numel(Ts2D) NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); ne00(:,:,it) = real(fftshift(ifft2((NE_),Nr,Nz))); ni00(:,:,it) = real(fftshift(ifft2((NI_),Nr,Nz))); phi (:,:,it) = real(fftshift(ifft2((PH_),Nr,Nz))); drphi(:,:,it) = real(fftshift(ifft2(1i*KR.*(PH_),Nr,Nz))); dr2phi(:,:,it)= real(fftshift(ifft2(-KR.^2.*(PH_),Nr,Nz))); dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz))); end % Building a version of phi only 5D sampling times PHI_Ts5D = zeros(Nkr,Nkz,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(Nkr,Nkz,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)),Nr,Nz))); end % Post processing disp('- post processing') E_pot = zeros(1,Ns2D); % Potential energy n^2 E_kin = zeros(1,Ns2D); % Kinetic energy grad(phi)^2 ExB = zeros(1,Ns2D); % ExB drift intensity \propto |\grad \phi| % gyrocenter and particle flux from real space GFlux_ri = zeros(1,Ns2D); % Gyrocenter flux Gamma = GFlux_zi = zeros(1,Ns2D); % Gyrocenter flux Gamma = GFlux_re = zeros(1,Ns2D); % Gyrocenter flux Gamma = GFlux_ze = zeros(1,Ns2D); % Gyrocenter flux Gamma = PFlux_ri = zeros(1,Ns5D); % Particle flux % gyrocenter and particle flux from fourier coefficients GFLUX_RI = real(squeeze(sum(sum(-1i*KZ.*Ni00.*conj(PHI),1),2)))*(2*pi/Nr/Nz)^2; PFLUX_RI = real(squeeze(sum(sum(-1i*KZ.*Np_i.*conj(PHI_Ts5D),1),2)))*(2*pi/Nr/Nz)^2; -phi_avg = zeros(1,Ns2D); % Time evol. of the norm of phi +phi_max = zeros(1,Ns2D); % Time evol. of the norm of phi Ne_norm = zeros(Npe,Nje,Ns5D); % Time evol. of the norm of Napj Ni_norm = zeros(Npi,Nji,Ns5D); % . Ddr = 1i*KR; Ddz = 1i*KZ; lapl = Ddr.^2 + Ddz.^2; for it = 1:numel(Ts2D) % Loop over 2D arrays NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); - phi_avg(it) = real(squeeze(PH_(ikr0,ikz0)))/2; + phi_max(it) = max(max(squeeze(phi(:,:,it)))); ExB(it) = max(max(max(abs(phi(3:end,:,it)-phi(1:end-2,:,it))/(2*dr))),max(max(abs(phi(:,3:end,it)-phi(:,1:end-2,it))'/(2*dz)))); GFlux_ri(it) = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz; GFlux_zi(it) = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz; GFlux_re(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz; GFlux_ze(it) = sum(sum(-ne00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz; 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)/Nkr/Nkz; Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkr/Nkz; % Particle flux PFlux_ri(it) = sum(sum(np_i(:,:,it).*dzphi(:,:,it2D)))*dr*dz/Lr/Lz; end %% Compute growth rate disp('- growth rate') % Find max value of transport (end of linear mode) [~,itmax] = max(GFlux_ri); tstart = 0.1 * Ts2D(itmax); tend = 0.9 * Ts2D(itmax); g_ = zeros(Nkr,Nkz); for ikr = 1:Nkr for ikz = 1:Nkz g_(ikr,ikz) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikr,ikz,:))),tstart,tend); end end [gmax,ikzmax] = max(g_(1,:)); kzmax = abs(kz(ikzmax)); Bohm_transport = ETAB/ETAN*gmax/kzmax^2; %% 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(221); + subplot(421); for ip = 1:Npe for ij = 1:Nje plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$N_e^{',num2str(Pe(ip)),num2str(Je(ij)),'}$']; 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(222) + subplot(423) for ip = 1:Npi for ij = 1:Nji plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$N_i^{',num2str(Pi(ip)),num2str(Ji(ij)),'}$']; 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}|$'); + grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); xlabel('$t c_s/R$') + subplot(222) + semilogy(Ts0D,GGAMMA_RI*(2*pi/Nr/Nz)^2); hold on; +% plot(Ts2D,GFLUX_RI) + plot(Ts0D,PGAMMA_RI*(2*pi/Nr/Nz)^2); +% plot(Ts5D,PFLUX_RI,'--'); + legend(['Gyro. flux';'Part. flux']); + grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma_{r,i}$') subplot(223) plot(kz,g_(1,:),'-','DisplayName','$\gamma$'); hold on; grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma R/c_s$'); %legend('show'); subplot(224) for ip = 1:Npi for ij = 1:Nji plt = @(x) squeeze(x(ip,ij,:)); plotname = '$\langle\phi\rangle_{r,z}(t)$'; clr = line_colors(min(ip,numel(line_colors(:,1))),:); lstyle = line_styles(min(ij,numel(line_styles))); - plot(Ts2D,phi_avg,'DisplayName',plotname,... + plot(Ts2D,phi_max,'DisplayName',plotname,... 'Color',clr,'LineStyle',lstyle{1}); hold on; end end - grid on; xlabel('$t c_s/R$'); ylabel('$\tilde\phi_{00}/2$'); %legend('show'); + grid on; xlabel('$t c_s/R$'); ylabel('$\max_{r,z}(\phi)$'); %legend('show'); % suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB)]); save_figure end if 0 %% Particle fluxes SCALING = Nkr*dkr * Nkz*dkz; fig = figure; FIGNAME = ['gamma',sprintf('_%.2d',JOBNUM),'_',PARAMS]; set(gcf, 'Position', [100, 100, 800, 300]) semilogy(Ts2D,GFLUX_RI, 'color', line_colors(2,:)); hold on plot(Ts5D,PFLUX_RI,'.', 'color', line_colors(2,:)); hold on plot(Ts2D,SCALING*GFlux_ri, 'color', line_colors(1,:)); hold on plot(Ts5D,SCALING*PFlux_ri,'.', 'color', line_colors(1,:)); hold on xlabel('$tc_{s0}/\rho_s$'); ylabel('$\Gamma_r$'); grid on title(['$\eta=',num2str(ETAB),'\quad',... '\nu_{',CONAME,'}=',num2str(NU),'$', ' $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']) legend('Gyro Flux','Particle flux', 'iFFT GFlux', 'iFFT PFlux')%'$\eta\gamma_{max}/k_{max}^2$') save_figure end if 0 %% Space time diagramm (fig 11 Ivanov 2020) fig = figure; FIGNAME = ['space_time_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600]) subplot(311) semilogy(Ts2D,GFLUX_RI,'-'); hold on plot(Ts5D,PFLUX_RI,'.'); hold on % plot(Ts2D,Bohm_transport*ones(size(Ts2D)),'--'); hold on ylabel('$\Gamma_r$'); grid on title(['$\eta=',num2str(ETAB),'\quad',... '\nu_{',CONAME,'}=',num2str(NU),'$']) legend(['$P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$'],'Particle flux')%'$\eta\gamma_{max}/k_{max}^2$') % set(gca,'xticklabel',[]) subplot(312) yyaxis left plot(Ts2D,squeeze(max(max((phi))))) ylabel('$\max \phi$') yyaxis right plot(Ts2D,squeeze(mean(max(dr2phi)))) ylabel('$s\sim\langle\partial_r^2\phi\rangle_z$'); grid on % set(gca,'xticklabel',[]) subplot(313) [TY,TX] = meshgrid(r,Ts2D); pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,:),2))'); set(pclr, 'edgecolor','none'); %colorbar; xlabel('$t c_s/R$'), ylabel('$r/\rho_s$') legend('$\langle\partial_r \phi\rangle_z$') save_figure end if 0 %% Photomaton : real space % FIELD = ni00; FNAME = 'ni'; % FIELD = ne00; FNAME = 'ne'; FIELD = phi; FNAME = 'phi'; tf = 200; [~,it1] = min(abs(Ts2D-tf)); tf = 600; [~,it2] = min(abs(Ts2D-tf)); tf =1000; [~,it3] = min(abs(Ts2D-tf)); tf =2000; [~,it4] = min(abs(Ts2D-tf)); fig = figure; FIGNAME = [FNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 400]) plt = @(x) x;%./max(max(x)); subplot(141) DATA = plt(FIELD(:,:,it1)); pclr = pcolor((RR),(ZZ),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$');set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it1))); subplot(142) DATA = plt(FIELD(:,:,it2)); pclr = pcolor((RR),(ZZ),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); subplot(143) DATA = plt(FIELD(:,:,it3)); pclr = pcolor((RR),(ZZ),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it3))); subplot(144) DATA = plt(FIELD(:,:,it4)); pclr = pcolor((RR),(ZZ),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it4))); % suptitle(['$\',FNAME,'$, $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),... % ', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']); save_figure end %% if 0 %% Show frame in kspace tf = 0; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf)); fig = figure; FIGNAME = ['krkz_',sprintf('t=%.0f',Ts2D(it2)),'_',PARAMS];set(gcf, 'Position', [100, 100, 700, 600]) subplot(221); plt = @(x) fftshift((abs(x)),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it2))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); legend('$|\hat\phi|$'); subplot(222); plt = @(x) fftshift(abs(x),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it2))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{00}|$'); subplot(223); plt = @(x) fftshift((abs(x)),2); FIELD = squeeze(Nipj(1,2,:,:,:)); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(FIELD(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{pj=01}|$'); subplot(224); plt = @(x) fftshift((abs(x)),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$');legend('$\hat S_i^{00}$'); save_figure end %% if 0 %% Ion moments max mode vs pj % tf = Ts2D(end-3); for tf = [] [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf)); % it2 = it2 + 1; fig = figure; FIGNAME = ['kmaxp_Nipj_',sprintf('t=%.2f',Ts2D(it2)),'_',PARAMS];set(gcf, 'Position', [100, 100, 700, 600]) plt = @(x) squeeze(max(abs(x),[],4)); % plt = @(x) squeeze(max(fftshift(abs(x),2),[],4)); for ij_ = 1:numel(Ji) subplot(100+numel(Ji)*10+ij_) pclr = imagesc(kr,Pi,plt(Nipj(:,ij_,:,:,it5))); xlabel('$k_r$'); if ij_ == 1 ylabel('$P$(max o. $k_z$)'); else yticks([]) end LEGEND = ['$|\hat n_i^{p',num2str(ij_-1),'}|$']; title(LEGEND); end save_figure end end %% t0 = 0; [~, it02D] = min(abs(Ts2D-t0)); [~, it05D] = min(abs(Ts5D-t0)); skip_ = 1; DELAY = 0.02*skip_; FRAMES_2D = it02D:skip_:numel(Ts2D); FRAMES_5D = it05D:skip_:numel(Ts5D); %% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if 0 %% Density ion GIFNAME = ['ni',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 1; FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D; FIELDNAME = '$n_i$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; create_gif end if 0 %% Phi real space GIFNAME = ['phi',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 1; FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D; FIELDNAME = '$\phi$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; create_gif end if 0 %% radial particle transport GIFNAME = ['gamma_r',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 1; FIELD = real(ni00.*dzphi); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D; FIELDNAME = '$\Gamma_r$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; create_gif end if 0 %% Phi fourier GIFNAME = ['FFT_phi',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 0; FIELD = ifftshift((abs(PHI)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D; FRAMES = FRAMES_2D; FIELDNAME = '$|\tilde\phi|$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$'; create_gif end if 0 %% phi @ z = 0 GIFNAME = ['phi_r0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; FIELD =(squeeze(real(phi(:,1,:)))); linestyle = '-.'; FRAMES = FRAMES_2D; X = (r); T = Ts2D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r); FIELDNAME = '$\phi(r=0)$'; XNAME = '$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(KR,2); Y = fftshift(KZ,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(KR,2); Y = fftshift(KZ,2); T = Ts2D; FIELDNAME = '$N_e^{00}$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$'; create_gif end if 0 %% kr vs P Si GIFNAME = ['Sip0_kr',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; plt = @(x) squeeze(max((abs(x)),[],4)); FIELD =plt(Sipj(:,1,:,:,:)); X = kr'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D; FIELDNAME = '$N_i^{p0}$'; XNAME = '$k_{max}\rho_s$'; YNAME = '$P$'; create_gif_imagesc end if 0 %% maxkz, kr vs p, for all Nipj over time GIFNAME = ['Nipj_kr',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; plt = @(x) squeeze(max((abs(x)),[],4)); FIELD = plt(Nipj); X = kr'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D; FIELDNAME = 'N_i'; XNAME = '$k_r\rho_s$'; YNAME = '$P$, ${k_z}^{max}$'; create_gif_5D end if 0 %% maxkr, kz vs p, for all Nipj over time GIFNAME = ['Nipj_kz',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; plt = @(x) fftshift(squeeze(max((abs(x)),[],3)),3); FIELD = plt(Nipj); X = sort(kz'); Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D; FIELDNAME = 'N_i'; XNAME = '$k_z\rho_s$'; YNAME = '$P$, ${k_r}^{max}$'; create_gif_5D end %% \ No newline at end of file diff --git a/wk/local_run.m b/wk/local_run.m index dfe9cce..9d897ea 100644 --- a/wk/local_run.m +++ b/wk/local_run.m @@ -1,55 +1,62 @@ addpath(genpath('../matlab')) % ... add %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 0.5; % Collision frequency ETAB = 0.6; % Magnetic gradient ETAN = 1.0; % Density gradient NU_HYP = 0.1; %% GRID PARAMETERS N = 50; % Frequency gridpoints (Nkr = N/2) L = 100; % Size of the squared frequency domain PMAXE = 10; % Highest electron Hermite polynomial degree JMAXE = 1; % Highest '' Laguerre '' PMAXI = 10; % Highest ion Hermite polynomial degree JMAXI = 1; % Highest '' Laguerre '' %% TIME PARAMETERS TMAX = 100; % Maximal time unit DT = 2e-2; % Time step -SPS0D = 1/DT; % Sampling per time unit for profiler +SPS0D = 1; % Sampling per time unit for profiler SPS2D = 1/2; % Sampling per time unit for 2D arrays SPS5D = 1/4; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints/10 RESTART = 0; % To restart from last checkpoint JOB2LOAD= 0; -%% OPTIONS +%% OPTIONS AND NAMING % SIMID = ['local_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation % SIMID = sprintf(SIMID,NU); % SIMID = 'test_init_phi'; % Name of the simulation SIMID = 'test_parallel_p'; % Name of the simulation 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) INIT_PHI= 1; % Start simulation with a noisy phi and moments +%% OUTPUTS +W_DOUBLE = 0; +W_GAMMA = 1; +W_PHI = 1; +W_NA00 = 1; +W_NAPJ = 1; +W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% unused KR0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. KREQ0 = 0; % put kr = 0 KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1) kmax = N*pi/L;% Highest fourier mode HD_CO = 0.5; % Hyper diffusivity cutoff ratio % kmaxcut = 2.5; MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient NOISE0 = 1.0e-5; TAU = 1.0; % e/i temperature ratio ETAT = 0.0; % Temperature gradient MU_P = 0.0/PMAXI^2; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_J = 0.0/JMAXI^3; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f %% Setup and file management setup system('rm fort.90'); \ No newline at end of file diff --git a/wk/marconi_run.m b/wk/marconi_run.m index 92e3531..48bab70 100644 --- a/wk/marconi_run.m +++ b/wk/marconi_run.m @@ -1,72 +1,79 @@ %clear all; addpath(genpath('../matlab')) % ... add %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% CLUSTER PARAMETERS CLUSTER.TIME = '01:00:00'; % allocation time hh:mm:ss CLUSTER.PART = 'dbg'; % dbg or prod CLUSTER.MEM = '16GB'; % Memory CLUSTER.JNAME = 'gamma_inf';% Job name NP_P = 2; % MPI processes along p NP_KR = 12; % MPI processes along kr %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 0.1; % Collision frequency ETAB = 0.6; % Magnetic gradient NU_HYP = 0.1; % Hyperdiffusivity coefficient %% GRID PARAMETERS -N = 050; % Frequency gridpoints (Nkr = N/2) -L = 100; % Size of the squared frequency domain -P = 10; % Electron and Ion highest Hermite polynomial degree -J = 01; % Electron and Ion highest Laguerre polynomial degree +N = 200; % Frequency gridpoints (Nkr = N/2) +L = 120; % Size of the squared frequency domain +P = 20; % Electron and Ion highest Hermite polynomial degree +J = 08; % 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 = 50; % Maximal time unit -DT = 2e-2; % Time step +TMAX = 200; % Maximal time unit +DT = 1e-2; % Time step SPS0D = 1; % Sampling per time unit for profiler -SPS2D = 1/2; % Sampling per time unit for 2D arrays -SPS5D = 1/10; % Sampling per time unit for 5D arrays +SPS2D = 1/10; % Sampling per time unit for 2D arrays +SPS5D = 0/50; % 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 = ['HeLaZ_v2.4_2_12_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation +SIMID = ['HeLaZ_v2.4_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation % SIMID = 'Marconi_test'; % 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) INIT_PHI= 1; % Start simulation with a noisy phi and moments +%% OUTPUTS +W_DOUBLE = 0; +W_GAMMA = 0; +W_PHI = 1; +W_NA00 = 1; +W_NAPJ = 0; +W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% fixed parameters (for current study) KR0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. KREQ0 = 0; % put kr = 0 KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 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_KR; Nnodes = ceil(Ntot/24); Nppn = Ntot/Nnodes; CLUSTER.NODES = num2str(Nnodes); % MPI process along p CLUSTER.NTPN = num2str(Nppn); % MPI process along kr CLUSTER.CPUPT = '1'; % CPU per task %% Run file management scripts setup write_sbash_marconi system('rm fort.90 setup_and_run.sh batch_script.sh'); disp('done');