diff --git a/.gitignore b/.gitignore index 201e2f7..963a466 100644 --- a/.gitignore +++ b/.gitignore @@ -1,40 +1,40 @@ .vimrc *.swp *~ *.asv *.m.orig *.DS_Store *.fig *.pdf *.eps *.mov *.mp4 *.emf *.pptx *.jpg *.jpeg *.mat *.dat *.mod *.h5 *.o *.a *.png logs/ results/ results_old/ mod/ obj/ bin/ -wk/fort.90 +wk/fort*.90 .directory checkpoint/ FM/ iCa/ *.out src/srcinfo.h src/srcinfo/srcinfo.h -fort.90 +fort*.90 local/ *.sh Gallery/ diff --git a/matlab/continue_run.m b/matlab/continue_run.m new file mode 100644 index 0000000..bf3c6a5 --- /dev/null +++ b/matlab/continue_run.m @@ -0,0 +1,92 @@ +%% Functions to modify preexisting fort.90 input file and launch on marconi +function [] = continue_run(outfilename) + EXECNAME = 'helaz_3.8'; + %% CLUSTER PARAMETERS + CLUSTER.PART = 'prod'; % dbg or prod + CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss + if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end; + CLUSTER.MEM = '64GB'; % Memory + CLUSTER.JNAME = 'HeLaZ';% Job name + NP_P = 2; % MPI processes along p + NP_KX = 24; % MPI processes along kx + % Compute processes distribution + Ntot = NP_P * NP_KX; + Nnodes = ceil(Ntot/48); + Nppn = Ntot/Nnodes; + CLUSTER.NODES = num2str(Nnodes); % MPI process along p + CLUSTER.NTPN = num2str(Nppn); % MPI process along kx + CLUSTER.CPUPT = '1'; % CPU per task + %% + RESDIR = ['../',outfilename(46:end-8),'/']; + BASIC.RESDIR = RESDIR; + FORT90 = [RESDIR,'fort.90']; + + % Read txt into cell A + fid = fopen(FORT90,'r'); + i = 1; + tline = fgetl(fid); + A{i} = tline; + while ischar(tline) + i = i+1; + tline = fgetl(fid); + A{i} = tline; + end + fclose(fid); + + % Find previous job2load + if( numel(A{5}) == numel(' RESTART = .false.') ) + A{5} = sprintf(' RESTART = .true.'); + J2L = 0; + else + line = A{39}; + line = line(end-2:end); + if(line(1) == '='); line = line(end); end; + J2L = str2num(line)+1; + end + % Change job 2 load in fort.90 + A{39} = [' job2load = ',num2str(J2L)]; + disp(A{39}) + % Change time step + line_= A{3}; + dt_old = str2num(line_(13:end)); + A{3} = [' dt = ',num2str(dt_old)]; + % Increase endtime + A{4} = [' tmax = 20000']; + % Change collision operator + line_= A{43}; + CO_old = str2num(line_(13:end)); + A{43} = [' CO = ',num2str(2)]; + % Put non linear term back + A{45} = [' NL_CLOS = -1']; + % change HD + line_= A{47}; + mu_old = str2num(line_(13:end)); + A{47} = [' mu = ',num2str(0)]; + % change L + line_= A{14}; + L_old = str2num(line_(12:end)); + A{14} = [' Lx = ',num2str(L_old)]; + A{16} = [' Ly = ',num2str(L_old)]; + % change eta N + line_= A{57}; + etan_old = str2num(line_(13:end)); + A{57} = [' eta_n = ',num2str(etan_old)]; + % change eta B + line_= A{59}; + etab_old = str2num(line_(13:end)); + A{59} = [' eta_B = ',num2str(etab_old)]; + % Rewrite fort.90 + fid = fopen('fort.90', 'w'); + for i = 1:numel(A) + if A{i+1} == -1 + fprintf(fid,'%s', A{i}); + break + else + fprintf(fid,'%s\n', A{i}); + end + end + % Copy fort.90 into marconi + write_sbash_marconi + % Launch the job + system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh'); +end diff --git a/matlab/load_marconi.m b/matlab/load_marconi.m index 05c0318..10406e8 100644 --- a/matlab/load_marconi.m +++ b/matlab/load_marconi.m @@ -1,25 +1,24 @@ function [ RESDIR ] = load_marconi( outfilename ) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here hostfolder = outfilename(1:end-8); hostfile = [hostfolder,'/out*']; MISCDIR = ['/misc/HeLaZ_outputs/',outfilename(46:end-8),'/']; RESDIR = ['../',outfilename(46:end-8),'/']; miscfolder = [MISCDIR,'.']; system(['mkdir -p ',miscfolder]); disp(['mkdir -p ',miscfolder]); resultfolder = [RESDIR,'.']; % SCP the output file from marconi to misc folder of SPCPC CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',miscfolder]; disp(CMD); system(CMD); % Load the fort.90 as well in misc folder - CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfolder,'/fort.90',' ',miscfolder]; + CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfolder,'/fort*.90',' ',miscfolder]; disp(CMD); system(CMD); % Put it also in the result directory - CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfolder,'/fort.90',' ',resultfolder]; + CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfolder,'/fort*.90',' ',resultfolder]; disp(CMD); system(CMD); end - diff --git a/matlab/setup.m b/matlab/setup.m index 168c965..6452cce 100644 --- a/matlab/setup.m +++ b/matlab/setup.m @@ -1,153 +1,153 @@ %% ________________________________________________________________________ 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.Nx = N; % x grid resolution GRID.Lx = L; % x length GRID.Ny = N * (1-KXEQ0) + KXEQ0; % y '' GRID.Ly = L * (1-KXEQ0); % y '' GRID.Nz = Nz; % z resolution GRID.q0 = q0; % q factor GRID.shear = shear; % shear GRID.eps = eps; % inverse aspect ratio % Model parameters MODEL.CO = CO; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) MODEL.CLOS = CLOS; MODEL.NL_CLOS = NL_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*KX0KH),'_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.INIT_ZF = INIT_ZF; INITIAL.wipe_turb = WIPE_TURB; INITIAL.wipe_zf = WIPE_ZF; if (INIT_BLOB); INITIAL.init_blob = '.true.'; else; INITIAL.init_blob = '.false.';end; INITIAL.init_background = (INIT_ZF>0)*ZF_AMP + BCKGD0; INITIAL.init_noiselvl = NOISE0; INITIAL.iseed = 42; INITIAL.mat_file = '''null'''; if (abs(CO) == 2) %Sugama operator INITIAL.mat_file = ['''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5''']; elseif (abs(CO) == 3) %pitch angle operator INITIAL.mat_file = ['''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5''']; elseif (CO == 4) % Full Coulomb GK % INITIAL.mat_file = ['''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5''']; INITIAL.mat_file = ['''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5''']; % INITIAL.mat_file = ['''../../../iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5''']; % INITIAL.mat_file = ['''../../../iCa/gk_coulomb_NFLR_6_P_4_J_2_N_75_kpm_6.0.h5''']; elseif (CO == -1) % DGDK disp('Warning, DGDK not debugged') end % Naming and creating input file switch abs(CO) case 0; CONAME = 'LB'; case 1; CONAME = 'DG'; case 2; CONAME = 'SG'; case 3; CONAME = 'PA'; case 4; CONAME = 'FC'; otherwise; CONAME ='UK'; end if (CO <= 0); CONAME = [CONAME,'DK']; else; CONAME = [CONAME,'GK']; 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,'_mu_%0.0e']; degngrad = sprintf(degngrad,[NU,MU]); if ~NON_LIN; degngrad = ['lin_',degngrad]; end if (Nz == 1) resolution = [num2str(GRID.Nx),'x',num2str(GRID.Ny/2),'_']; gridname = ['L_',num2str(L),'_']; else resolution = [num2str(GRID.Nx),'x',num2str(GRID.Ny/2),'x',num2str(GRID.Nz),'_']; gridname = ['L_',num2str(L),'_q0_',num2str(q0),'_']; end if (exist('PREFIX','var') == 0); PREFIX = []; end; if (exist('SUFFIX','var') == 0); SUFFIX = []; end; PARAMS = [PREFIX,resolution,gridname,degngrad,SUFFIX]; BASIC.RESDIR = [SIMDIR,PARAMS,'/']; BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',SIMDIR(4:end),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_3d = floor(1.0/SPS3D/DT); OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT); OUTPUTS.nsave_cp = floor(1.0/SPSCP/DT); 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_HF; OUTPUTS.write_hf = '.true.'; else; OUTPUTS.write_hf = '.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; if W_DENS; OUTPUTS.write_dens = '.true.'; else; OUTPUTS.write_dens = '.false.';end; if W_TEMP; OUTPUTS.write_temp = '.true.'; else; OUTPUTS.write_temp = '.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 if ~exist(BASIC.MISCDIR, 'dir') mkdir(BASIC.MISCDIR) 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 +if JOB2LOAD>=0 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 d23fdec..b4bf4f0 100644 --- a/matlab/write_fort90.m +++ b/matlab/write_fort90.m @@ -1,89 +1,88 @@ function [INPUT] = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC) % Write the input script "fort.90" with desired parameters -INPUT = 'fort.90'; +INPUT = ['fort_',sprintf('%2.2d',OUTPUTS.job2load+1),'.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,[' 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,[' Nx = ', num2str(GRID.Nx),'\n']); fprintf(fid,[' Lx = ', num2str(GRID.Lx),'\n']); fprintf(fid,[' Ny = ', num2str(GRID.Ny),'\n']); fprintf(fid,[' Ly = ', num2str(GRID.Ly),'\n']); fprintf(fid,[' Nz = ', num2str(GRID.Nz),'\n']); fprintf(fid,[' q0 = ', num2str(GRID.q0),'\n']); fprintf(fid,[' shear = ', num2str(GRID.shear),'\n']); fprintf(fid,[' eps = ', num2str(GRID.eps),'\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_3d = ', num2str(OUTPUTS.nsave_3d),'\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_hf = ', OUTPUTS.write_hf,'\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,[' write_dens = ', OUTPUTS.write_dens,'\n']); fprintf(fid,[' write_temp = ', OUTPUTS.write_temp,'\n']); -fprintf(fid,[' resfile0 = ', OUTPUTS.resfile0,'\n']); -fprintf(fid,[' rstfile0 = ', OUTPUTS.rstfile0,'\n']); -fprintf(fid,[' job2load = ', num2str(OUTPUTS.job2load),'\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,[' NL_CLOS = ', num2str(MODEL.NL_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,[' INIT_ZF = ', num2str(INITIAL.INIT_ZF),'\n']); fprintf(fid,[' WIPE_ZF = ', num2str(INITIAL.wipe_zf),'\n']); fprintf(fid,[' WIPE_TURB = ', num2str(INITIAL.wipe_turb),'\n']); fprintf(fid,[' INIT_BLOB = ', INITIAL.init_blob,'\n']); fprintf(fid,[' init_background = ', num2str(INITIAL.init_background),'\n']); fprintf(fid,[' init_noiselvl = ', num2str(INITIAL.init_noiselvl),'\n']); fprintf(fid,[' iseed = ', num2str(INITIAL.iseed),'\n']); fprintf(fid,[' mat_file = ', INITIAL.mat_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,'/.']); +system(['cp fort*.90 ',BASIC.RESDIR,'/.']); end diff --git a/matlab/write_sbash_daint.m b/matlab/write_sbash_daint.m index 62a48ac..536ff2f 100644 --- a/matlab/write_sbash_daint.m +++ b/matlab/write_sbash_daint.m @@ -1,57 +1,57 @@ % Write the input script "fort.90" with desired parameters INPUT = 'setup_and_run.sh'; fid = fopen(INPUT,'wt'); fprintf(fid,[... '#!/bin/bash\n',... 'mkdir -p $SCRATCH/HeLaZ/wk\n',... ... 'cd $SCRATCH/HeLaZ/wk/\n',... ... 'mkdir -p ', BASIC.RESDIR,'\n',... 'cd ',BASIC.RESDIR,'\n',... -'cp $HOME/HeLaZ/wk/fort.90 .\n',... +'cp $HOME/HeLaZ/wk/fort*.90 .\n',... 'cp $HOME/HeLaZ/wk/batch_script.sh .\n',... ... 'jid=$(sbatch batch_script.sh)\n',... 'echo $jid\n',... 'echo to check output log :\n',... 'echo tail -f $SCRATCH/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS,'/out.txt']); fclose(fid); system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']); % Write the sbatch script INPUT = 'batch_script.sh'; fid = fopen(INPUT,'wt'); fprintf(fid,[... '#!/bin/bash\n',... '#SBATCH --job-name="',CLUSTER.JNAME,'"\n',... '#SBATCH --time=', CLUSTER.TIME,'\n',... '#SBATCH --nodes=', CLUSTER.NODES,'\n',... '#SBATCH --cpus-per-task=', CLUSTER.CPUPT,'\n',... '#SBATCH --ntasks-per-node=', CLUSTER.NTPN,'\n',... '#SBATCH --ntasks-per-core=', CLUSTER.NTPC,'\n',... '#SBATCH --mem=', CLUSTER.MEM,'\n',... '#SBATCH --error=err.txt\n',... '#SBATCH --output=out.txt\n',... '#SBATCH --account="s882"\n',... '#SBATCH --constraint=mc\n',... '#SBATCH --hint=nomultithread\n',... '#SBATCH --partition=',CLUSTER.PART,'\n',... 'export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK\n',... ...% '#SBATCH --job-name=',PARAMS,'\n\n',... 'module purge\n',... 'module load PrgEnv-intel\n',... 'module load cray-hdf5-parallel\n',... 'module load cray-mpich\n',... 'module load craype-x86-skylake\n',... 'module load cray-fftw\n',... 'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KX)]); %'srun ./../../../bin/helaz']); fclose(fid); system(['cp batch_script.sh ',BASIC.RESDIR,'/.']); -system('scp {fort.90,setup_and_run.sh,batch_script.sh} ahoffman@ela.cscs.ch:HeLaZ/wk'); \ No newline at end of file +system('scp {fort*.90,setup_and_run.sh,batch_script.sh} ahoffman@ela.cscs.ch:HeLaZ/wk'); diff --git a/matlab/write_sbash_marconi.m b/matlab/write_sbash_marconi.m index d4fde78..461b0c5 100644 --- a/matlab/write_sbash_marconi.m +++ b/matlab/write_sbash_marconi.m @@ -1,45 +1,45 @@ % Write the input script "fort.90" with desired parameters INPUT = 'setup_and_run.sh'; fid = fopen(INPUT,'wt'); fprintf(fid,[... '#!/bin/bash\n',... 'mkdir -p $CINECA_SCRATCH/HeLaZ/wk\n',... ... 'cd $CINECA_SCRATCH/HeLaZ/wk/\n',... ... 'mkdir -p ', BASIC.RESDIR,'\n',... 'cd ',BASIC.RESDIR,'\n',... -'cp $HOME/HeLaZ/wk/fort.90 .\n',... +'cp $HOME/HeLaZ/wk/fort*.90 .\n',... 'cp $HOME/HeLaZ/wk/batch_script.sh .\n',... ... -'sbatch batch_script.sh\n',... -'echo tail -f $CINECA_SCRATCH/HeLaZ',BASIC.RESDIR(3:end),'out.txt']); +SBATCH_CMD,... +'echo tail -f $CINECA_SCRATCH/HeLaZ',BASIC.RESDIR(3:end),'out']); fclose(fid); system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']); % Write the sbatch script INPUT = 'batch_script.sh'; fid = fopen(INPUT,'wt'); fprintf(fid,[... '#!/bin/bash\n',... '#SBATCH --job-name=',CLUSTER.JNAME,'\n',... '#SBATCH --time=', CLUSTER.TIME,'\n',... '#SBATCH --nodes=', CLUSTER.NODES,'\n',... '#SBATCH --cpus-per-task=', CLUSTER.CPUPT,'\n',... '#SBATCH --ntasks-per-node=', CLUSTER.NTPN,'\n',... '#SBATCH --mem=', CLUSTER.MEM,'\n',... -'#SBATCH --error=err.txt\n',... -'#SBATCH --output=out.txt\n',... +'#SBATCH --error=err',num2str(JOB2LOAD+1),'.txt\n',... +'#SBATCH --output=out_',num2str(JOB2LOAD+1),'.txt\n',... '#SBATCH --account=FUA35_TSVVT421\n',... '#SBATCH --partition=skl_fua_',CLUSTER.PART,'\n',... 'module load autoload hdf5 fftw\n',... -'srun --cpu-bind=cores ./../../../bin/',EXECNAME,' ',num2str(NP_P),' ',num2str(NP_KX)]); +'srun --cpu-bind=cores ./../../../bin/',EXECNAME,' ',num2str(NP_P),' ',num2str(NP_KX),' ',num2str(JOB2LOAD+1)]); fclose(fid); system(['cp batch_script.sh ',BASIC.RESDIR,'/.']); -system('scp {fort.90,setup_and_run.sh,batch_script.sh} ahoffman@login.marconi.cineca.it:/marconi/home/userexternal/ahoffman/HeLaZ/wk > trash.txt'); +system('scp {fort*.90,setup_and_run.sh,batch_script.sh} ahoffman@login.marconi.cineca.it:/marconi/home/userexternal/ahoffman/HeLaZ/wk > trash.txt'); system('rm trash.txt'); \ No newline at end of file diff --git a/src/tesend.F90 b/src/tesend.F90 index 16ef076..7e356e9 100644 --- a/src/tesend.F90 +++ b/src/tesend.F90 @@ -1,70 +1,70 @@ SUBROUTINE tesend ! Test for run completion USE basic use prec_const IMPLICIT NONE LOGICAL :: mlend, mlexist REAL :: tnow INTEGER :: ncheck_stop = 100 CHARACTER(len=*), PARAMETER :: stop_file = 'mystop' !________________________________________________________________________________ ! 1. Some processors had set nlend CALL mpi_allreduce(nlend, mlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr) IF( mlend ) THEN nlend = .TRUE. crashed = .TRUE. IF (my_id .EQ. 0) WRITE(*,'(/a)') 'rhs are NaN/Inf' IF (my_id .EQ. 0) WRITE(*,*) 'Run terminated at cstep=',cstep RETURN END IF !________________________________________________________________________________ ! 2. Test on NRUN nlend = step .GT. nrun IF ( nlend ) THEN WRITE(*,'(/a)') 'NRUN steps done' RETURN END IF !________________________________________________________________________________ ! 3. Test on TMAX nlend = time .GT. tmax IF ( nlend ) THEN IF (my_id .EQ. 0) WRITE(*,'(/a)') 'TMAX reached' RETURN END IF ! ! !________________________________________________________________________________ ! 4. Test on run time CALL cpu_time(tnow) mlend = (1.2*(tnow-start)) .GT. maxruntime CALL mpi_allreduce(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr) IF ( nlend ) THEN IF(my_id.EQ.0) WRITE(*,'(/a)') 'Max run time reached' RETURN END IF !________________________________________________________________________________ - ! 5. NRUN modified throught "stop file" + ! 5. NRUN modified through "stop file" ! IF( (my_id .EQ. 0) .AND. (MOD(cstep, ncheck_stop) == 0) ) THEN INQUIRE(file=stop_file, exist=mlexist) IF( mlexist ) THEN OPEN(lu_stop, file=stop_file) mlend = mlexist ! Send stop status asa the file exists WRITE(*,'(/a,i4,a)') 'Stop file found -> finishing..' CLOSE(lu_stop, status='delete') END IF END IF CALL mpi_allreduce(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr) ! RETURN ! END SUBROUTINE tesend diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m index 0bdc0cf..dc2ba46 100644 --- a/wk/analysis_3D.m +++ b/wk/analysis_3D.m @@ -1,310 +1,310 @@ addpath(genpath('../matlab')) % ... add addpath(genpath('../matlab/plots')) % ... add outfile =''; %% Directory of the simulation -if 0% Local results +if 1% Local results outfile =''; outfile =''; % outfile ='artificial_ZF_freeze/sim_A'; % outfile ='simulation_B/cw_FCGK_kp_3.0'; -% outfile ='nonlin_FCGK/150x75_L_200_P_4_J_2_eta_0.6_nu_1e-01_FCGK_mu_0e+00'; +outfile ='nonlin_FCGK/150x75_L_200_P_4_J_2_eta_0.6_nu_1e-01_FCGK_mu_0e+00'; % outfile ='nonlin_PAGK/100x50_L_200_P_4_J_2_eta_0.6_nu_1e-01_PAGK_mu_0e+00'; % outfile ='nonlin_FCGK/100x50_L_200_P_4_J_2_eta_0.6_nu_1e-01_FCGK_mu_0e+00'; -outfile ='simulation_A'; +% outfile ='simulation_A'; % outfile ='simulation_B/cw_SGGK_like_species'; % outfile ='simulation_A/CO_damping_SGGK'; % outfile ='simulation_A/cw_DGGK_eta_0.5'; % outfile ='simulation_B/cw_DGGK'; BASIC.RESDIR = ['../results/',outfile,'/']; BASIC.MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/']; system(['mkdir -p ',BASIC.MISCDIR]); CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD); system(CMD); else% Marconi results outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt'; % outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt'; BASIC.RESDIR = ['../',outfile(46:end-8),'/']; BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/']; end %% Load the results % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 20; % JOBNUMMIN = 07; JOBNUMMAX = 20; % For CO damping sim A compile_results %Compile the results from first output found to JOBNUMMAX if existing %% Post-processing post_processing %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 0 %% Time evolutions and growth rate plot_time_evolution_and_gr end if 1 %% Space time diagramm (fig 11 Ivanov 2020) TAVG_0 = 1000; TAVG_1 = 2000; % Averaging times duration plot_radial_transport_and_shear end if 0 %% Space time diagramms cmax = 0.00001 % max of the colorbar for transport tstart = 0; tend = Ts3D(end); % time window plot_space_time_diagrams end if 0 %% |phi_k|^2 spectra (Kobayashi 2015 fig 3) % tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window tstart = 2500; tend = 2500; % tstart = 10000; tend = 12000; % Chose the field to plot % FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}'; % FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}' FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi'; % FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:76,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x'; LOGSCALE = 1; TRENDS = 1; NORMALIZED = 0; plot_kperp_spectrum end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options t0 =000; iz = 1; ix = 1; iy = 1; skip_ =4; DELAY = 2e-3*skip_; [~, it03D] = min(abs(Ts3D-t0)); FRAMES_3D = it03D:skip_:numel(Ts3D); [~, it05D] = min(abs(Ts5D-t0)); FRAMES_5D = it05D:skip_:numel(Ts5D); T = Ts3D; FRAMES = FRAMES_3D; INTERP = 0; NORMALIZED = 0; CONST_CMAP = 0;% Gif options % Field to plot % FIELD = dens_e; NAME = 'ne'; FIELDNAME = 'n_e'; % FIELD = dens_i; NAME = 'ni'; FIELDNAME = 'n_i'; % FIELD = dens_e-Z_n_e; NAME = 'ne_NZ';FIELDNAME = 'n_e^{NZ}'; FIELD = dens_i-Z_n_i; NAME = 'ni_NZ';FIELDNAME = 'n_i^{NZ}'; % FIELD = temp_e; NAME = 'Te'; FIELDNAME = 'n_i'; % FIELD = temp_i; NAME = 'Ti'; FIELDNAME = 'n_i'; % FIELD = temp_e-Z_T_e; NAME = 'Te_NZ';FIELDNAME = 'T_e^{NZ}'; % FIELD = temp_i-Z_T_i; NAME = 'Ti_NZ';FIELDNAME = 'T_i^{NZ}'; % FIELD = ne00; NAME = 'ne00'; FIELDNAME = 'n_e^{00}'; % FIELD = ni00; NAME = 'ni00'; FIELDNAME = 'n_i^{00}'; % FIELD = phi; NAME = 'phi'; FIELDNAME = '\phi'; % FIELD = Z_phi; NAME = 'Zphi'; FIELDNAME = '\phi_Z'; % FIELD = Gamma_x; NAME = 'Gamma_x'; FIELDNAME = '\Gamma_x'; % Sliced % plt = @(x) real(x(ix, :, :,:)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; % plt = @(x) real(x( :,iy, :,:)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; % plt = @(x) real(x( :, :,iz,:)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; % Averaged % plt = @(x) mean(x,1); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; % plt = @(x) mean(x,2); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; plt = @(x) mean(x,3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELD = squeeze(plt(FIELD)); % Naming GIFNAME = [NAME,sprintf('_%.2d',JOBNUM),'_',PARAMS]; % Create movie (gif or mp4) create_gif % create_mov end if 0 %% Photomaton : real space % Chose the field to plot % FIELD = ni00; FNAME = 'ni00'; FIELDLTX = 'n_i^{00}'; % FIELD = ne00; FNAME = 'ne00'; FIELDLTX = 'n_e^{00}' % FIELD = dens_i; FNAME = 'ni'; FIELDLTX = 'n_i'; % FIELD = dens_e; FNAME = 'ne'; FIELDLTX = 'n_e'; % FIELD = dens_e-Z_n_e; FNAME = 'ne_NZ'; FIELDLTX = 'n_e^{NZ}'; % FIELD = dens_i-Z_n_i; FNAME = 'ni_NZ'; FIELDLTX = 'n_i^{NZ}'; % FIELD = temp_i; FNAME = 'Ti'; FIELDLTX = 'T_i'; % FIELD = temp_e; FNAME = 'Te'; FIELDLTX = 'T_e'; % FIELD = phi; FNAME = 'phi'; FIELDLTX = '\phi'; % FIELD = Z_phi-phi; FNAME = 'phi_NZ'; FIELDLTX = '\phi^{NZ}'; FIELD = Z_phi; FNAME = 'phi_Z'; FIELDLTX = '\phi^{Z}'; % FIELD = Gamma_x; FNAME = 'Gamma_x'; FIELDLTX = '\Gamma_x'; % FIELD = dens_e-Z_n_e-(Z_phi-phi); FNAME = 'Non_adiab_part'; FIELDLTX = 'n_e^{NZ}-\phi^{NZ}'; % Chose when to plot it tf = [9800 10000]; % Sliced ix = 1; iy = 1; iz = 1; % plt = @(x,it) real(x(ix, :, :,it)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(x=',num2str(round(x(ix))),')'] % plt = @(x,it) real(x( :,iy, :,it)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(y=',num2str(round(y(iy))),')'] plt = @(x,it) real(x( :, :,iz,it)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = [FIELDLTX,'(z=',num2str(round(z(iz)/pi)),'\pi)'] % Averaged % plt = @(x,it) mean(x(:,:,:,it),1); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; FIELDLTX = ['\langle ',FIELDLTX,'\rangle_x'] % plt = @(x,it) mean(x(:,:,:,it),2); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; FIELDLTX = ['\langle ',FIELDLTX,'\rangle_y'] % plt = @(x,it) mean(x(:,:,:,it),3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; FIELDLTX = ['\langle ',FIELDLTX,'\rangle_z'] % TNAME = []; fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 350]) plt_2 = @(x) x;%./max(max(x)); for i_ = 1:numel(tf) [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))]; subplot(1,numel(tf),i_) DATA = plt_2(squeeze(plt(FIELD,it))); pclr = pcolor((X),(Y),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) colormap(bluewhitered); %caxis([-30,30]); xlabel(latexize(XNAME)); ylabel(latexize(YNAME));set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts3D(it))); end legend(latexize(FIELDLTX)); save_figure end if 0 %% Photomaton : k space % Chose the field to plot % FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}'; % FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}' FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi'; % FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x'; % FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e'; % Chose when to plot it tf = 1000:500:3000; % tf = 8000; % Sliced ix = 1; iy = 1; iz = 1; % plt = @(x,it) abs(x(ix, :, :,it)); X = KY_YZ; Y = KZ_YZ; XNAME = 'k_y'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(k_x=',num2str(round(kx(ix))),')']; % plt = @(x,it) abs(x( :,iy, :,it)); X = KX_XZ; Y = KZ_XZ; XNAME = 'k_x'; YNAME = 'z'; FIELDLTX = [FIELDLTX,'(k_y=',num2str(round(ky(iy))),')']; plt = @(x,it) abs(x( :, :,iz,it)); X = KX_XY; Y = KY_XY; XNAME = 'k_x'; YNAME = 'k_y'; FIELDLTX = [FIELDLTX,'(z=',num2str((z(iz)/pi)),'\pi)']; % TNAME = []; fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 300*numel(tf), 500]) plt_2 = @(x) (fftshift(x,2)); for i_ = 1:numel(tf) [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))]; subplot(1,numel(tf),i_) DATA = plt_2(squeeze(plt(FIELD,it))); pclr = pcolor(fftshift(X,2),fftshift(Y,2),DATA); set(pclr, 'edgecolor','none');pbaspect([0.5 1 1]) caxis([0 1]*5e3); xlabel(latexize(XNAME)); ylabel(latexize(YNAME)); if(i_ > 1); set(gca,'ytick',[]); end; title(sprintf('$t c_s/R=%.0f$',Ts3D(it))); end legend(latexize(FIELDLTX)); save_figure end if 0 %% TAVG_0 = 1000; TAVG_1 = 5000; % Averaging times duration ZF_fourier_analysis end if 0 %% plot_param_evol end if 0 %% Radial shear profile tf = 3000+[900:20:1100]; ymax = 0; figure for i_ = 1:numel(tf) [~,it] = min(abs(Ts3D-tf(i_))); data = squeeze((mean(dxphi(:,:,1,it),2))); plot(x,data,'Displayname',sprintf('$t c_s/R=%.0f$',Ts3D(it))); hold on; ymax = max([ymax abs(min(data)) abs(max(data))]); end xlim([min(x), max(x)]); ylim(1.2*[-1 1]*ymax); xlabel('$x/\rho_s$'); ylabel('$v_{E\times B,x}$'); grid on end if 1 %% zonal vs nonzonal energies for phi(t) trange = 200:Ns3D; Ephi_Z = zeros(1,Ns3D); Ephi_NZ_kgt0 = zeros(1,Ns3D); Ephi_NZ_kgt1 = zeros(1,Ns3D); Ephi_NZ_kgt2 = zeros(1,Ns3D); high_k_phi = zeros(1,Ns3D); for it = 1:numel(Ts3D) % Ephi_NZ(it) = sum(sum(((KY~=0).*abs(PHI(:,:,1,it)).^2))); % Ephi_Z(it) = sum(sum(((KY==0).*abs(PHI(:,:,1,it)).^2))); [amp,ikzf] = max(abs((kx~=0).*PHI(:,1,1,it))); % Ephi_NZ(it) = sum(sum(((KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2))); Ephi_NZ_kgt0(it) = sum(sum(((sqrt(KX.^2+KY.^2)>0.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2))); Ephi_NZ_kgt1(it) = sum(sum(((sqrt(KX.^2+KY.^2)>1.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2))); Ephi_NZ_kgt2(it) = sum(sum(((sqrt(KX.^2+KY.^2)>2.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2))); % Ephi_Z(it) = kx(ikzf)^2*abs(PHI(ikzf,1,1,it)).^2; Ephi_Z(it) = sum(sum(((KX~=0).*(KY==0).*(KX.^2).*abs(PHI(:,:,1,it)).^2))); % Ephi_NZ(it) = sum(sum(((KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)))-Ephi_Z(it); high_k_phi(it) = abs(PHI(18,18,1,it)).^2; % kperp sqrt(2) % high_k_phi(it) = abs(PHI(40,40,1,it)).^2;% kperp 3.5 end pltx = @(x) x-x(1); plty = @(x) x./max(squeeze(x)); fig = figure; FIGNAME = ['ZF_turb_energy_vs_time_',PARAMS]; set(gcf, 'Position', [100, 100, 1400, 500]) subplot(121) % yyaxis left semilogy(pltx(Ts3D(trange)),plty(Ephi_Z(trange)),'DisplayName',['Zonal, ',CONAME]); hold on; % yyaxis right semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt0(trange)),'DisplayName',['NZ, $k_p>0$, ',CONAME]); semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt1(trange)),'DisplayName',['NZ, $k_p>1$, ',CONAME]); semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt2(trange)),'DisplayName',['NZ, $k_p>2$, ',CONAME]); % semilogy(pltx(Ts0D),plty(PGAMMA_RI),'DisplayName',['$\Gamma_x$, ',CONAME]); title('Energy'); legend('show') xlabel('$t c_s/R$'); grid on;% xlim([0 500]); subplot(122) plot(plty(Ephi_Z(trange)),plty(Ephi_NZ_kgt0(trange))); title('Phase space'); legend(CONAME) xlabel('$E_Z$'); ylabel('$E_{NZ}$'); grid on;% xlim([0 500]); % % subplot(133) % % semilogy(pltx(Ts0D),plty(abs(PGAMMA_RI)*SCALE)); % for ik = [10 20 30] % semilogy(pltx(Ts3D(trange)),plty(abs(squeeze(PHI(ik,ik,1,trange))).^2),'DisplayName',[CONAME,', kp=',num2str(sqrt(kx(ik)^2+ky(ik)^2))]); hold on % end % title('High kperp damping'); legend('show'); % xlabel('$t c_s/R$'); grid on;% xlim([0 500]); end if 0 %% Conservation laws Nxmax = Nx/2; Nymax = Ny/2; mflux_x_i = squeeze(sum((Gamma_x( 1,1:Nxmax,:)+Gamma_x( 1,2:Nxmax+1,:))/2,2)./sum(Gamma_x( 1,1:Nxmax,:))); mflux_x_o = squeeze(sum((Gamma_x( Nxmax,1:Nxmax,:)+Gamma_x( Nxmax,2:Nxmax+1,:))/2,2)./sum(Gamma_x( Nxmax,1:Nxmax,:))); mflux_y_i = squeeze(sum((Gamma_y(1:Nxmax, 1,:)+Gamma_y(2:Nxmax+1, 1,:))/2,1)./sum(Gamma_y(1:Nxmax, 1,:))); mflux_y_o = squeeze(sum((Gamma_y(1:Nxmax, Nymax,:)+Gamma_y(2:Nxmax+1, Nymax,:))/2,1)./sum(Gamma_y(1:Nxmax, Nymax,:))); mass_cons = mflux_x_i - mflux_x_o + mflux_y_i - mflux_y_o; %% figure plt = @(x) squeeze(mean(mean(x(:,:,1,:),1),2)); subplot(211) plot(Ts3D,plt(dens_e+dens_i),'DisplayName','$\delta n_e + \delta n_i$'); hold on; plot(Ts3D,plt(ne00+ni00),'DisplayName','$\delta n_e^{00} + \delta n_i^{00}$'); hold on; plot(Ts3D,plt(temp_e+temp_i),'DisplayName','$\delta T_e + \delta T_i$'); hold on; legend('show'); grid on; xlim([Ts3D(1) Ts3D(end)]) subplot(212); plot(Ts3D,mass_cons*(2*pi/Nx/Ny)^2,'DisplayName','in-out'); hold on % plot(Ts3D,squeeze(mflux_x_i),'DisplayName','Flux i x'); % plot(Ts3D,squeeze(mflux_x_o),'DisplayName','Flux o x'); % plot(Ts3D,squeeze(mflux_y_i),'DisplayName','Flux i y'); % plot(Ts3D,squeeze(mflux_y_o),'DisplayName','Flux o y'); legend('show'); grid on; xlim([Ts3D(1) Ts3D(end)]); %ylim([-0.1, 2]*mean(mflux_x_i)) end \ No newline at end of file diff --git a/wk/local_run.m b/wk/local_run.m index 7de7329..06bcd40 100644 --- a/wk/local_run.m +++ b/wk/local_run.m @@ -1,78 +1,74 @@ addpath(genpath('../matlab')) % ... add %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 0.1; % Collision frequency ETAN = 1.0/0.6; % Density gradient drive (R/Ln) NU_HYP = 0.0; %% GRID AND GEOMETRY PARAMETERS N = 150; % Frequency gridpoints (Nkx = N/2) L = 200; % Size of the squared frequency domain Nz = 1; % number of perpendicular planes (parallel grid) q0 = 1.0; % safety factor shear = 0.0; % magnetic shear eps = 0.0; % inverse aspect ratio P = 4; J = 2; %% TIME PARAMETERS TMAX = 100; % Maximal time unit DT = 2e-2; % Time step SPS0D = 1; % Sampling per time unit for profiler SPS2D = 1; % Sampling per time unit for 2D arrays SPS3D = 1; % Sampling per time unit for 3D arrays SPS5D = 1/20; % 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; +JOB2LOAD= -1; %% OPTIONS AND NAMING % Collision operator % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; 4 : Coulomb; +/- for GK/DK) -CO = 4; +CO = 1; CLOS = 0; % Closure model (0: =0 truncation) NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS) -SIMID = 'nonlin_FCGK'; % Name of the simulation -% SIMID = 'test_3D'; % Name of the simulation +% SIMID = 'nonlin_FCGK'; % Name of the simulation +SIMID = 'test'; % Name of the simulation % SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation NON_LIN = 1; % activate non-linearity (is cancelled if KXEQ0 = 1) % INIT options INIT_ZF = 0; ZF_AMP = 0.0; INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0; %% OUTPUTS -W_DOUBLE = 0; -W_GAMMA = 1; -W_PHI = 1; -W_NA00 = 1; -W_NAPJ = 1; -W_SAPJ = 0; -W_DENS = 1; -W_TEMP = 1; +W_DOUBLE = 1; +W_GAMMA = 1; W_HF = 1; +W_PHI = 1; W_NA00 = 1; +W_DENS = 1; W_TEMP = 1; +W_NAPJ = 1; W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% unused PMAXE = P; % Highest electron Hermite polynomial degree JMAXE = J; % Highest '' Laguerre '' PMAXI = P; % Highest ion Hermite polynomial degree JMAXI = J; % Highest '' Laguerre '' KERN = 0; % Kernel model (0 : GK) KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. KXEQ0 = 0; % put kx = 0 KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; 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 +MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient NOISE0 = 1.0e-5; +BCKGD0 = 0.0; % Init background TAU = 1.0; % e/i temperature ratio ETAT = 0.0; % Temperature gradient ETAB = 1.0; % Magnetic gradient (1.0 to set R=LB) INIT_PHI= 1; % Start simulation with a noisy phi and moments MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f %% Setup and file management setup -system('rm fort.90'); +system('rm fort*.90'); outfile = [BASIC.RESDIR,'out.txt']; disp(outfile); diff --git a/wk/marconi_run.m b/wk/marconi_run.m index 577e6b5..9c2b2b6 100644 --- a/wk/marconi_run.m +++ b/wk/marconi_run.m @@ -1,107 +1,118 @@ clear all; addpath(genpath('../matlab')) % ... add SUBMIT = 1; % To submit the job automatically +CHAIN = 1; % To chain jobs (CHAIN = n will launch n jobs in chain) % EXECNAME = 'helaz_dbg'; - EXECNAME = 'helaz_3.8'; + EXECNAME = 'helaz_3.81'; for ETAN = [1/0.6] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% CLUSTER PARAMETERS CLUSTER.PART = 'prod'; % dbg or prod % CLUSTER.PART = 'dbg'; -CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss +CLUSTER.TIME = '20:00:00'; % allocation time hh:mm:ss if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end; CLUSTER.MEM = '128GB'; % Memory CLUSTER.JNAME = 'HeLaZ';% Job name NP_P = 2; % MPI processes along p -NP_KX = 24; % MPI processes along kx +NP_KX = 48; % MPI processes along kx %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS -NU = 0.5; % Collision frequency +NU = 0.1; % Collision frequency ETAN = 1.0/0.6; % Density gradient drive (R/Ln) NU_HYP = 0.0; %% GRID PARAMETERS N = 300; % Frequency gridpoints (Nkx = N/2) L = 120; % Size of the squared frequency domain Nz = 1; % number of perpendicular planes (parallel grid) q0 = 1.0; % q factor () shear = 0.0; % magnetic shear eps = 0.0; % inverse aspect ratio P = 8; J = 4; %% TIME PARAMETERS TMAX = 10000; % Maximal time unit -DT = 5e-3; % Time step +DT = 8e-3; % Time step SPS0D = 1; % Sampling per time unit for profiler SPS2D = 1; % Sampling per time unit for 2D arrays SPS3D = 1/2; % Sampling per time unit for 3D arrays SPS5D = 1/100; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints/10 -RESTART = 1; % To restart from last checkpoint -JOB2LOAD= 0; +JOB2LOAD= -1; % start from t=0 if <0, else restart from outputs_$job2load %% OPTIONS AND NAMING % Collision operator % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK) -CO = 2; +CO = 3; CLOS = 0; % Closure model (0: =0 truncation) NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS) -% SIMID = 'test_3D_marconi'; % Name of the simulation -SIMID = 'simulation_B'; % Name of the simulation +% SIMID = 'test_chained_job'; % Name of the simulation +SIMID = 'simulation_A'; % Name of the simulation % SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation NON_LIN = 1; % activate non-linearity (is cancelled if KXEQ0 = 1) % INIT options INIT_ZF = 0; ZF_AMP = 0.0; INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0; %% OUTPUTS W_DOUBLE = 1; -W_GAMMA = 1; -W_PHI = 1; -W_NA00 = 1; -W_NAPJ = 1; -W_SAPJ = 0; -W_DENS = 1; -W_TEMP = 1; +W_GAMMA = 1; W_HF = 1; +W_PHI = 1; W_NA00 = 1; +W_DENS = 1; W_TEMP = 1; +W_NAPJ = 1; W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% unused PMAXE = P; % Highest electron Hermite polynomial degree JMAXE = J; % Highest '' Laguerre '' PMAXI = P; % Highest ion Hermite polynomial degree JMAXI = J; % Highest '' Laguerre '' KERN = 0; % Kernel model (0 : GK) KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. KXEQ0 = 0; % put kx = 0 KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; 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; BCKGD0 = 0.0; % Init background TAU = 1.0; % e/i temperature ratio ETAT = 0.0; % Temperature gradient ETAB = 1.0; % Magnetic gradient (1.0 to set R=LB) INIT_PHI= 1; % Start simulation with a noisy phi and moments MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f % Compute processes distribution Ntot = NP_P * NP_KX; Nnodes = ceil(Ntot/48); Nppn = Ntot/Nnodes; CLUSTER.NODES = num2str(Nnodes); % MPI process along p CLUSTER.NTPN = num2str(Nppn); % MPI process along kx CLUSTER.CPUPT = '1'; % CPU per task %% Run file management scripts setup +SBATCH_CMD = 'sbatch batch_script.sh\n'; write_sbash_marconi -system('rm fort.90 setup_and_run.sh batch_script.sh'); +system('rm fort*.90 setup_and_run.sh batch_script.sh'); if(mod(NP_P*NP_KX,48)~= 0) disp('WARNING : unused cores (ntot cores must be a 48 multiple)'); end if(SUBMIT) - system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh'); + [~,job_info_] = system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh'); + disp(job_info_); + jobid_ = job_info_(21:27); + if(CHAIN>0) + for CHAIN_IDX = 1:CHAIN + SBATCH_CMD = ['sbatch --dependency=afterok:',jobid_,' batch_script.sh\n']; + disp(SBATCH_CMD); + JOB2LOAD= JOB2LOAD+1; + setup + write_sbash_marconi + [~,job_info_] = system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh'); + jobid_ = job_info_(21:27); + end + end end disp('done'); end