diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m index f10777e..686686f 100644 --- a/matlab/plot/plot_cosol_mat.m +++ b/matlab/plot/plot_cosol_mat.m @@ -1,180 +1,181 @@ % MAT = results.iCa; % figure % suptitle('FCGK,P=6,J=3'); % subplot(221) % imagesc(log(abs(MAT))); % title('log abs') % subplot(222) % imagesc(imag(MAT)); colormap(bluewhitered) % title('imag'); colorbar % subplot(223) % imagesc(imag(MAT)>0); % title('imag$>$0'); % subplot(224) % imagesc(imag(MAT)<0); % title('imag$<$0'); % % %% SGGK P_ = 20; J_ = 10; mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'; SGGK_self = h5read(mat_file_name,'/00000/Caapj/Ciipj'); SGGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT'); SGGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT'); figure MAT = 1i*SGGK_self; suptitle('SGGK ii,P=20,J=10, k=0'); % MAT = 1i*SGGK_ei; suptitle('SGGK ei,P=20,J=10'); % MAT = 1i*SGGK_ie; suptitle('SGGK ie,P=20,J=10'); subplot(221) imagesc(abs(MAT)); title('log abs') subplot(222) imagesc(imag(MAT)); colormap(bluewhitered) title('imag'); colorbar subplot(223) imagesc(imag(MAT)>0); title('imag$>$0'); subplot(224) imagesc(imag(MAT)<0); title('imag$<$0'); %% PAGK P_ = 20; J_ = 10; mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'; PAGK_self = h5read(mat_file_name,'/00000/Caapj/Ceepj'); % PAGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT'); % PAGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT'); figure MAT = 1i*PAGK_self; suptitle('PAGK ii,P=20,J=10, k=0'); subplot(221) imagesc(abs(MAT)); title('log abs') subplot(222) imagesc(imag(MAT)); colormap(bluewhitered) title('imag'); colorbar subplot(223) imagesc(imag(MAT)>0); title('imag$>$0'); subplot(224) imagesc(imag(MAT)<0); title('imag$<$0'); %% FCGK P_ = 4; J_ = 2; -% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5'; mat_file_name = '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5'; +% mat_file_name = '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5'; % mat_file_name = '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5'; -kp = 2.0; +kp = 0.0; kp_a = h5read(mat_file_name,'/coordkperp'); [~,matidx] = min(abs(kp_a-kp)); matidx = sprintf('%5.5i',matidx); FCGK_self = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']); FCGK_ei = h5read(mat_file_name,['/',matidx,'/Ceipj/CeipjF'])+h5read(mat_file_name,'/00000/Ceipj/CeipjT'); FCGK_ie = h5read(mat_file_name,['/',matidx,'/Ciepj/CiepjF'])+h5read(mat_file_name,'/00000/Ciepj/CiepjT'); figure MAT = 1i*FCGK_self; suptitle(['FCGK ii,P=6,J=3, k=',num2str(kp),' (',matidx,')']); % MAT = 1i*FCGK_ei; suptitle('FCGK ei,P=20,J=10'); % MAT = 1i*FCGK_ie; suptitle('FCGK ie,P=20,J=10'); subplot(221) imagesc(abs(MAT)); title('log abs') subplot(222) imagesc(imag(MAT)); colormap(bluewhitered) title('imag'); colorbar subplot(223) imagesc(imag(MAT)>0); title('imag$>$0'); subplot(224) imagesc(imag(MAT)<0); title('imag$<$0'); %% Eigenvalue spectrum analysis if 0 %% mfns = {... '/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',... % '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',... -% '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_4.h5',... -% '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5',... -% '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5',... - '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_30.h5',... +% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_4.h5',... +% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12.h5',... +% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5',... +% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5',... + '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5',... % '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',... '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',... % '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5',... % '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5',... }; CONAME_A = {... 'SG 20 10',... % 'PA 20 10',... % 'FC 10 5 NFLR 4',... % 'FC 10 5 NFLR 12',... % 'FC 10 5 NFLR 12 k<2', ... - 'FC 10 5 NFLR 30', ... + 'LD 6 3 NFLR 20', ... % 'FC 4 2 NFLR 6',... 'FC 4 2 NFLR 12', ... % 'Hacked SG A',... % 'Hacked SG B',... }; figure for j_ = 1:numel(mfns) mat_file_name = mfns{j_}; disp(mat_file_name); kp_a = h5read(mat_file_name,'/coordkperp'); % kp_a = kp_a(kp_a<=3); gmax = zeros(size(kp_a)); wmax = zeros(size(kp_a)); for idx_ = 0:numel(kp_a)-1 matidx = sprintf('%5.5i',idx_); MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']); gmax(idx_+1) = max(abs(real(eig(MAT)))); wmax(idx_+1) = max(abs(imag(eig(MAT)))); end subplot(121) plot(kp_a,gmax,'DisplayName',CONAME_A{j_}); hold on; subplot(122) plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on; end subplot(121) legend('show'); grid on; ylim([0,100]); xlabel('$k_\perp$'); ylabel('$\gamma_{max}$ from Eig(iCa)') subplot(122) legend('show'); grid on; xlabel('$k_\perp$'); ylabel('$\omega_{max}$ from Eig(iCa)') end %% Van Kampen plot if 0 %% kperp= 1.5; -mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',... - '/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',... - '/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',... - '/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',... - '/home/ahoffman/HeLaZ/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5',... +mfns = {'/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',... + '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',... + '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',... + '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',... + '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5',... }; -CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'FC 10 5 NFLR 12'}; +CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'LD 6 3 NFLR 12'}; grow = {}; puls = {}; for j_ = 1:numel(mfns) mat_file_name = mfns{j_}; disp(mat_file_name); kp_a = h5read(mat_file_name,'/coordkperp'); [~,idx_] = min(abs(kp_a-kperp)); matidx = sprintf('%5.5i',idx_); MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']); grow{j_} = real(eig(MAT)); puls{j_} = imag(eig(MAT)); end figure for j_ = 1:numel(mfns) % plot(puls{j_}, grow{j_},'o','DisplayName',CONAME_A{j_}); hold on; plot(grow{j_},'o','DisplayName',CONAME_A{j_}); hold on; end legend('show'); grid on; title(['$k_\perp=$',num2str(kperp)]); xlabel('$\omega$ from Eig(iCa)'); ylabel('$\gamma$ from Eig(iCa)') end \ No newline at end of file diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m index 7b95c93..d548c76 100644 --- a/matlab/plot/plot_metric.m +++ b/matlab/plot/plot_metric.m @@ -1,30 +1,30 @@ function [ fig ] = plot_metric( data ) names = {'Jacobian','gradxB','gradyB','gradzB','gradz_coeff',... 'gxx','gxy','gyy','gyz','gzz','hatB','hatR','hatZ'}; geo_arrays = zeros(2,data.Nz,numel(names)); for i_ = 1:numel(names) namae = names{i_}; geo_arrays(:,:,i_) = h5read(data.outfilenames{end},['/data/metric/',namae])'; end fig = figure; subplot(311) for i = 1:5 plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on; end - xlim([min(data.z),max(data.z)]); legend('show'); title('MoNoLiT geometry'); + xlim([min(data.z),max(data.z)]); legend('show'); title('GYACOMO geometry'); subplot(312) for i = 6:10 plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on; end xlim([min(data.z),max(data.z)]); legend('show'); subplot(313) for i = 11:13 plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on; end xlim([min(data.z),max(data.z)]); legend('show'); end diff --git a/matlab/plot/statistical_transport_averaging.m b/matlab/plot/statistical_transport_averaging.m index eb284c6..3a7cc92 100644 --- a/matlab/plot/statistical_transport_averaging.m +++ b/matlab/plot/statistical_transport_averaging.m @@ -1,45 +1,50 @@ function [ fig ] = statistical_transport_averaging( data, options ) scale = data.scale; Trange = options.T; [~,it0] = min(abs(Trange(1)-data.Ts0D)); [~,it1] = min(abs(Trange(end)-data.Ts0D)); gamma = data.PGAMMA_RI(it0:it1)*scale; Qx = data.HFLUX_X(it0:it1)*scale; dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1; % if ~dt_const % disp('DT not const on given interval'); % else Ntot = (it1-it0)+1; transp_seg_avg = 1:Ntot; transp_seg_std = 1:Ntot; for Np = 1:Ntot % Loop on the number of segments transp_seg_avg(Np) = mean(gamma(1:Np)); transp_seg_std(Np) = std(gamma(1:Np)); end time_seg = (data.Ts0D(it0:it1)-data.Ts0D(it0)); - + + fig = 0; +if options.NPLOTS > 0 fig = figure; subplot(211) plot(time_seg,transp_seg_avg,'-'); hold on; xlabel('Averaging time'); ylabel('$\langle\Gamma_x\rangle_{\tau}$'); legend(['$\Gamma_x^\infty=$',sprintf('%2.2e',transp_seg_avg(end))]) title(sprintf('Transport averaging from t=%2.2f',data.Ts0D(it0))); for Np = 1:Ntot % Loop on the number of segments transp_seg_avg(Np) = mean(Qx(1:Np)); end subplot(212) plot(time_seg,transp_seg_avg,'-'); hold on; xlabel('Averaging time'); ylabel('$\langle Q_x\rangle_{\tau}$'); legend(['$Q_x^\infty=$',sprintf('%2.2e',transp_seg_avg(end))]) - +end +Gx_infty_avg = mean(gamma); +Gx_infty_std = std (gamma); +disp(['G_x=',sprintf('%2.2e',Gx_infty_avg),'+-',sprintf('%2.2e',Gx_infty_std)]); end diff --git a/matlab/setup.m b/matlab/setup.m index 966b6fb..2661841 100644 --- a/matlab/setup.m +++ b/matlab/setup.m @@ -1,186 +1,187 @@ %% _______________________________________________________________________ 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 = NX; % x grid resolution GRID.Lx = LX; % x length GRID.Nexc = NEXC; % to extend Lx when s>0 GRID.Ny = NY; % y '' GRID.Ly = LY; % y '' GRID.Nz = NZ; % z resolution GRID.Npol = NPOL; % z resolution if SG; GRID.SG = '.true.'; else; GRID.SG = '.false.';end; % Geometry GEOM.geom = ['''',GEOMETRY,'''']; GEOM.q0 = Q0; % q factor GEOM.shear = SHEAR; % shear GEOM.eps = EPS; % inverse aspect ratio % Model parameters MODEL.CLOS = CLOS; MODEL.NL_CLOS = NL_CLOS; MODEL.LINEARITY = ['''',LINEARITY,'''']; MODEL.KIN_E = KIN_E; if KIN_E; MODEL.KIN_E = '.true.'; else; MODEL.KIN_E = '.false.';end; MODEL.beta = BETA; MODEL.mu_x = MU_X; MODEL.mu_y = MU_Y; MODEL.N_HD = N_HD; MODEL.mu_z = MU_Z; 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 = SIGMA_E; 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.K_Ni = K_Ni; MODEL.K_Ne = K_Ne; MODEL.K_Ti = K_Ti; MODEL.K_Te = K_Te; MODEL.GradB = GRADB; % Magnetic gradient MODEL.CurvB = CURVB; % Magnetic curvature MODEL.lambdaD = LAMBDAD; % Collision parameters COLL.collision_model = ['''',CO,'''']; if (GKCO); COLL.gyrokin_CO = '.true.'; else; COLL.gyrokin_CO = '.false.';end; if (ABCO); COLL.interspecies = '.true.'; else; COLL.interspecies = '.false.';end; COLL.mat_file = '''null'''; switch CO case 'SG' COLL.mat_file = '''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'''; % COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5'''; % COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5'''; case 'LR' COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'''; case 'LD' % COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5'''; -% COLL.mat_file = '''../../../iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5'''; - COLL.mat_file = '''../../../iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_30.h5'''; +% COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5'''; +% COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5'''; + COLL.mat_file = '''../../../iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5'''; end COLL.coll_kcut = COLL_KCUT; % Time integration and intialization parameters TIME_INTEGRATION.numerical_scheme = '''RK4'''; INITIAL.INIT_OPT = ['''',INIT_OPT,'''']; INITIAL.ACT_ON_MODES = ['''',ACT_ON_MODES,'''']; INITIAL.init_background = BCKGD0; INITIAL.init_noiselvl = NOISE0; INITIAL.iseed = 42; % Naming and creating input file CONAME = CO; if GKCO CONAME = [CONAME,'GK']; else CONAME = [CONAME,'DK']; end if ~ABCO CONAME = [CONAME,'aa']; end if (CLOS == 0); CLOSNAME = 'Trunc.'; elseif(CLOS == 1); CLOSNAME = 'Clos. 1'; elseif(CLOS == 2); CLOSNAME = 'Clos. 2'; end % Hermite-Laguerre degrees naming if (PMAXE == PMAXI) && (JMAXE == JMAXI) HLdeg_ = ['_',num2str(PMAXE+1),'x',num2str(JMAXE+1)]; else HLdeg_ = ['_Pe_',num2str(PMAXE+1),'_Je_',num2str(JMAXE+1),... '_Pi_',num2str(PMAXI+1),'_Ji_',num2str(JMAXI+1)]; end % temp. dens. drives drives_ = []; if abs(K_Ni) > 0; drives_ = [drives_,'_kN_',num2str(K_Ni)]; end; if abs(K_Ti) > 0; drives_ = [drives_,'_kT_',num2str(K_Ti)]; end; % collision coll_ = ['_nu_%1.1e_',CONAME]; coll_ = sprintf(coll_,NU); % nonlinear lin_ = []; if ~LINEARITY; lin_ = '_lin'; end adiabe_ = []; if ~KIN_E; adiabe_ = '_adiabe'; end % resolution and boxsize res_ = [num2str(GRID.Nx),'x',num2str(GRID.Ny)]; if (LX ~= LY) geo_ = ['_Lx_',num2str(LX),'_Ly_',num2str(LY)]; else geo_ = ['_L_',num2str(LX)]; end if (NZ > 1) %3D case res_ = [res_,'x',num2str(NZ)]; if abs(Q0) > 0 geo_ = [geo_,'_q0_',num2str(Q0)]; end if abs(EPS) > 0 geo_ = [geo_,'_e_',num2str(EPS)]; end if abs(SHEAR) > 0 geo_ = [geo_,'_s_',num2str(SHEAR)]; end end switch GEOMETRY case 'circular' geo_ = [geo_,'_circ_']; end % put everything together in the param character chain u_ = '_'; % underscore variable PARAMS = [res_,HLdeg_,geo_,drives_,coll_,lin_,adiabe_]; 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 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); 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.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,GEOM,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC); nproc = 1; MAKE = 'cd ..; make; cd wk'; % system(MAKE); %% disp(['Set up ',SIMID]); disp([res_,geo_,HLdeg_]); if JOB2LOAD>=0 disp(['- restarting from JOBNUM = ',num2str(JOB2LOAD)]); else disp(['- starting from T = 0']); end diff --git a/src/diagnose.F90 b/src/diagnose.F90 index bfd9e9b..0763216 100644 --- a/src/diagnose.F90 +++ b/src/diagnose.F90 @@ -1,646 +1,648 @@ SUBROUTINE diagnose(kstep) ! Diagnostics, writing simulation state to disk USE basic USE diagnostics_par IMPLICIT NONE INTEGER, INTENT(in) :: kstep CALL cpu_time(t0_diag) ! Measuring time !! Basic diagnose loop for reading input file, displaying advancement and ending IF ((kstep .EQ. 0)) THEN INQUIRE(unit=lu_in, name=input_fname) CLOSE(lu_in) ENDIF IF (kstep .GE. 0) THEN ! Terminal info IF (MOD(cstep, INT(1.0/dt)) == 0 .AND. (my_id .EQ. 0)) THEN WRITE(*,"(F6.0,A,F6.0)") time,"/",tmax ENDIF ELSEIF (kstep .EQ. -1) THEN CALL cpu_time(finish) ! Display computational time cost IF (my_id .EQ. 0) CALL display_h_min_s(finish-start) END IF !! Specific diagnostic calls CALL diagnose_full(kstep) CALL cpu_time(t1_diag); tc_diag = tc_diag + (t1_diag - t0_diag) END SUBROUTINE diagnose SUBROUTINE init_outfile(comm,file0,file,fid) USE diagnostics_par, ONLY : write_doubleprecision, diag_par_outputinputs, input_fname USE basic, ONLY : my_id, jobnum, lu_in, basic_outputinputs USE grid, ONLY : grid_outputinputs USE geometry, ONLY : geometry_outputinputs USE model, ONLY : model_outputinputs USE collision, ONLY : coll_outputinputs USE initial_par, ONLY : initial_outputinputs USE time_integration,ONLY : time_integration_outputinputs USE futils, ONLY : creatf, creatg, creatd, attach, putfile IMPLICIT NONE !input INTEGER, INTENT(IN) :: comm CHARACTER(len=256), INTENT(IN) :: file0 CHARACTER(len=256), INTENT(OUT) :: file INTEGER, INTENT(OUT) :: fid CHARACTER(len=256) :: str,fname INCLUDE 'srcinfo.h' ! Writing output filename WRITE(file,'(a,a1,i2.2,a3)') TRIM(file0) ,'_',jobnum,'.h5' ! 1.1 Initial run ! Main output file creation IF (write_doubleprecision) THEN CALL creatf(file, fid, real_prec='d', mpicomm=comm) ELSE CALL creatf(file, fid, mpicomm=comm) END IF IF (my_id .EQ. 0) WRITE(*,'(3x,a,a)') TRIM(file), ' created' ! basic data group CALL creatg(fid, "/data", "data") ! File group CALL creatg(fid, "/files", "files") CALL attach(fid, "/files", "jobnum", jobnum) ! Add the code info and parameters to the file WRITE(str,'(a,i2.2)') "/data/input" CALL creatd(fid, 0,(/0/),TRIM(str),'Input parameters') CALL attach(fid, TRIM(str), "version", VERSION) !defined in srcinfo.h CALL attach(fid, TRIM(str), "branch", BRANCH) !defined in srcinfo.h CALL attach(fid, TRIM(str), "author", AUTHOR) !defined in srcinfo.h CALL attach(fid, TRIM(str), "execdate", EXECDATE) !defined in srcinfo.h CALL attach(fid, TRIM(str), "host", HOST) !defined in srcinfo.h CALL basic_outputinputs(fid,str) CALL grid_outputinputs(fid, str) CALL geometry_outputinputs(fid, str) CALL diag_par_outputinputs(fid, str) CALL model_outputinputs(fid, str) CALL coll_outputinputs(fid, str) CALL initial_outputinputs(fid, str) CALL time_integration_outputinputs(fid, 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 CALL putfile(fid, TRIM(str), TRIM(input_fname),ionode=0) END SUBROUTINE init_outfile SUBROUTINE diagnose_full(kstep) USE basic USE grid USE diagnostics_par USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf, putarrnd USE array USE model USE initial_par USE fields USE time_integration USE parallel USE prec_const USE collision, ONLY: coll_outputinputs USE geometry IMPLICIT NONE INTEGER, INTENT(in) :: kstep INTEGER, parameter :: BUFSIZE = 2 INTEGER :: rank = 0 INTEGER :: dims(1) = (/0/) INTEGER :: cp_counter = 0 CHARACTER(len=256) :: str,test_ !____________________________________________________________________________ ! 1. Initial diagnostics IF ((kstep .EQ. 0)) THEN CALL init_outfile(comm0, resfile0,resfile,fidres) ! 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_clos", "cumulative closure computation time") CALL creatd(fidres, 0, dims, "/profiler/Tc_ghost", "cumulative communication time") CALL creatd(fidres, 0, dims, "/profiler/Tc_coll", "cumulative collision computation 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_checkfield", "cumulative checkfield computation time") CALL creatd(fidres, 0, dims, "/profiler/Tc_diag", "cumulative sym computation time") CALL creatd(fidres, 0, dims, "/profiler/Tc_process", "cumulative process 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") ! Grid info CALL creatg(fidres, "/data/grid", "Grid data") CALL putarr(fidres, "/data/grid/coordkx", kxarray_full, "kx*rho_s0", ionode=0) CALL putarr(fidres, "/data/grid/coordky", kyarray_full, "ky*rho_s0", ionode=0) CALL putarr(fidres, "/data/grid/coordz", zarray_full, "z/R", ionode=0) CALL putarr(fidres, "/data/grid/coordp_e" , parray_e_full, "p_e", ionode=0) CALL putarr(fidres, "/data/grid/coordj_e" , jarray_e_full, "j_e", ionode=0) CALL putarr(fidres, "/data/grid/coordp_i" , parray_i_full, "p_i", ionode=0) CALL putarr(fidres, "/data/grid/coordj_i" , jarray_i_full, "j_i", ionode=0) ! Metric info CALL creatg(fidres, "/data/metric", "Metric data") CALL putarrnd(fidres, "/data/metric/gxx", gxx(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/gxy", gxy(izs:ize,0:1), (/1, 1, 1/)) + CALL putarrnd(fidres, "/data/metric/gxz", gxz(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/gyy", gyy(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/gyz", gyz(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/gzz", gzz(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/hatR", hatR(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/hatZ", hatZ(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/hatB", hatB(izs:ize,0:1), (/1, 1, 1/)) + CALL putarrnd(fidres, "/data/metric/hatB_NL", hatB_NL(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/gradxB", gradxB(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/gradyB", gradyB(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/gradzB", gradzB(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/Jacobian", Jacobian(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff(izs:ize,0:1), (/1, 1, 1/)) CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/1, 1, 3/)) CALL putarrnd(fidres, "/data/metric/kernel_i", kernel_i(ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/ 1, 2, 4/)) ! 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") IF(KIN_E) THEN CALL creatd(fidres, rank, dims, "/data/var0d/gflux_re", "Radial gyro electron transport") CALL creatd(fidres, rank, dims, "/data/var0d/pflux_re", "Radial part electron transport") ENDIF ENDIF IF (write_hf) THEN CALL creatd(fidres, rank, dims, "/data/var0d/hflux_xi", "Radial part ion heat flux") IF(KIN_E) THEN CALL creatd(fidres, rank, dims, "/data/var0d/hflux_xe", "Radial part electron heat flux") ENDIF ENDIF IF (cstep==0) THEN iframe0d=0 ENDIF CALL attach(fidres,"/data/var0d/" , "frames", iframe0d) END IF ! var2d group (??) 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 (cstep==0) THEN iframe2d=0 ENDIF CALL attach(fidres,"/data/var2d/" , "frames", iframe2d) END IF ! var3d group (electro. pot., Ni00 moment) IF (nsave_3d .GT. 0) THEN CALL creatg(fidres, "/data/var3d", "3d profiles") CALL creatd(fidres, rank, dims, "/data/var3d/time", "Time t*c_s/R") CALL creatd(fidres, rank, dims, "/data/var3d/cstep", "iteration number") IF (write_phi) CALL creatg(fidres, "/data/var3d/phi", "phi") IF (write_phi) CALL creatg(fidres, "/data/var3d/psi", "psi") IF (write_Na00) THEN IF(KIN_E)& CALL creatg(fidres, "/data/var3d/Ne00", "Ne00") CALL creatg(fidres, "/data/var3d/Ni00", "Ni00") IF(KIN_E)& CALL creatg(fidres, "/data/var3d/Nepjz", "Nepjz") CALL creatg(fidres, "/data/var3d/Nipjz", "Nipjz") ENDIF IF (write_dens) THEN IF(KIN_E)& CALL creatg(fidres, "/data/var3d/dens_e", "dens_e") CALL creatg(fidres, "/data/var3d/dens_i", "dens_i") ENDIF IF (write_fvel) THEN IF(KIN_E) THEN CALL creatg(fidres, "/data/var3d/upar_e", "upar_e") CALL creatg(fidres, "/data/var3d/uper_e", "uper_e") ENDIF CALL creatg(fidres, "/data/var3d/upar_i", "upar_i") CALL creatg(fidres, "/data/var3d/uper_i", "uper_i") ENDIF IF (write_temp) THEN IF(KIN_E) THEN CALL creatg(fidres, "/data/var3d/Tper_e", "Tper_e") CALL creatg(fidres, "/data/var3d/Tpar_e", "Tpar_e") CALL creatg(fidres, "/data/var3d/temp_e", "temp_e") ENDIF CALL creatg(fidres, "/data/var3d/Tper_i", "Tper_i") CALL creatg(fidres, "/data/var3d/Tpar_i", "Tpar_i") CALL creatg(fidres, "/data/var3d/temp_i", "temp_i") ENDIF IF (cstep==0) THEN iframe3d=0 ENDIF CALL attach(fidres,"/data/var3d/" , "frames", iframe3d) END IF ! var5d group (moments) IF (nsave_5d .GT. 0) THEN 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 IF(KIN_E)& CALL creatg(fidres, "/data/var5d/moments_e", "moments_e") CALL creatg(fidres, "/data/var5d/moments_i", "moments_i") ENDIF IF (write_Sapj) THEN IF(KIN_E)& CALL creatg(fidres, "/data/var5d/Sepj", "Sepj") CALL creatg(fidres, "/data/var5d/Sipj", "Sipj") ENDIF IF (cstep==0) THEN iframe5d=0 END IF CALL attach(fidres,"/data/var5d/" , "frames", iframe5d) END IF ENDIF !_____________________________________________________________________________ ! 2. Periodic diagnostics ! IF (kstep .GE. 0) THEN ! 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 ! empty in our case ! 2.3 3d profiles IF (nsave_3d .GT. 0) THEN IF (MOD(cstep, nsave_3d) == 0) THEN CALL diagnose_3d ! Looks at the folder if the file check_phi exists and spits a snapshot ! of the current electrostatic potential in a basic text file CALL spit_snapshot_check ENDIF ENDIF ! 2.4 5d 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 attach(fidres, "/data/input","cpu_time",finish-start) ! make a checkpoint at last timestep if not crashed IF(.NOT. crashed) THEN IF(my_id .EQ. 0) write(*,*) 'Saving last state' IF (nsave_5d .GT. 0) & CALL diagnose_5d ENDIF ! Close all diagnostic files CALL mpi_barrier(MPI_COMM_WORLD, ierr) CALL closef(fidres) END IF END SUBROUTINE diagnose_full !!-------------- Auxiliary routines -----------------!! SUBROUTINE diagnose_0d USE basic USE futils, ONLY: append, attach, getatt USE diagnostics_par USE prec_const USE processing USE model, ONLY: KIN_E 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_clos", tc_clos,ionode=0) CALL append(fidres, "/profiler/Tc_ghost", tc_ghost,ionode=0) CALL append(fidres, "/profiler/Tc_coll", tc_coll,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_checkfield",tc_checkfield,ionode=0) CALL append(fidres, "/profiler/Tc_diag", tc_diag,ionode=0) CALL append(fidres, "/profiler/Tc_process", tc_process,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) IF(KIN_E) THEN CALL compute_radial_electron_transport CALL append(fidres, "/data/var0d/gflux_re",gflux_re,ionode=0) CALL append(fidres, "/data/var0d/pflux_re",pflux_re,ionode=0) ENDIF ENDIF IF (write_hf) THEN CALL compute_radial_ion_heatflux CALL append(fidres, "/data/var0d/hflux_xi",hflux_xi,ionode=0) IF(KIN_E) THEN CALL compute_radial_electron_heatflux CALL append(fidres, "/data/var0d/hflux_xe",hflux_xe,ionode=0) ENDIF ENDIF END SUBROUTINE diagnose_0d SUBROUTINE diagnose_3d USE basic USE futils, ONLY: append, getatt, attach, putarrnd, putarr USE fields USE array USE grid, ONLY: ikxs,ikxe, ikys,ikye, Nkx, Nky, local_nkx, ikx, iky, ips_e, ips_i USE time_integration USE diagnostics_par USE prec_const USE processing USE model, ONLY: KIN_E IMPLICIT NONE INTEGER :: i_, root, world_rank, world_size CHARACTER(256) :: dset_name CALL append(fidres, "/data/var3d/time", time,ionode=0) CALL append(fidres, "/data/var3d/cstep", real(cstep,dp),ionode=0) CALL getatt(fidres, "/data/var3d/", "frames",iframe3d) iframe3d=iframe3d+1 CALL attach(fidres,"/data/var3d/" , "frames", iframe3d) IF (write_phi) CALL write_field3d_kykxz(phi (ikys:ikye,ikxs:ikxe,izs:ize), 'phi') IF (write_phi) CALL write_field3d_kykxz(psi (ikys:ikye,ikxs:ikxe,izs:ize), 'psi') IF (write_Na00) THEN IF(KIN_E)THEN IF (CONTAINS_ip0_e) & Ne00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_e(ip0_e,ij0_e,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel) CALL write_field3d_kykxz(Ne00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ne00') ENDIF IF (CONTAINS_ip0_i) & Ni00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_i(ip0_i,ij0_i,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel) CALL write_field3d_kykxz(Ni00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ni00') ! CALL compute_Napjz_spectrum ! IF(KIN_E) & ! CALL write_field3d_pjz_e(Nepjz(ips_e:ipe_e,ijs_e:ije_e,izs:ize), 'Nepjz') ! CALL write_field3d_pjz_i(Nipjz(ips_i:ipe_i,ijs_i:ije_i,izs:ize), 'Nipjz') ENDIF !! Fuid moments IF (write_dens .OR. write_fvel .OR. write_temp) & CALL compute_fluid_moments IF (write_dens) THEN IF(KIN_E)& CALL write_field3d_kykxz(dens_e(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_e') CALL write_field3d_kykxz(dens_i(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_i') ENDIF IF (write_fvel) THEN IF(KIN_E)& CALL write_field3d_kykxz(upar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_e') CALL write_field3d_kykxz(upar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_i') IF(KIN_E)& CALL write_field3d_kykxz(uper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_e') CALL write_field3d_kykxz(uper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_i') ENDIF IF (write_temp) THEN IF(KIN_E)& CALL write_field3d_kykxz(Tpar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_e') CALL write_field3d_kykxz(Tpar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_i') IF(KIN_E)& CALL write_field3d_kykxz(Tper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_e') CALL write_field3d_kykxz(Tper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_i') IF(KIN_E)& CALL write_field3d_kykxz(temp_e(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_e') CALL write_field3d_kykxz(temp_i(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_i') ENDIF CONTAINS SUBROUTINE write_field3d_kykxz(field, text) USE parallel, ONLY : gather_xyz IMPLICIT NONE COMPLEX(dp), DIMENSION(ikys:ikye,ikxs:ikxe, izs:ize), INTENT(IN) :: field CHARACTER(*), INTENT(IN) :: text COMPLEX(dp), DIMENSION(1:Nky,1:Nkx,1:Nz) :: field_full CHARACTER(256) :: dset_name WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d IF (num_procs .EQ. 1) THEN ! no data distribution CALL putarr(fidres, dset_name, field(ikys:ikye,ikxs:ikxe, izs:ize), ionode=0) ELSEIF(GATHERV_OUTPUT) THEN ! output using one node (gatherv) CALL gather_xyz(field(ikys:ikye,1:Nkx,izs:ize),field_full(1:Nky,1:Nkx,1:Nz)) CALL putarr(fidres, dset_name, field_full(1:Nky,1:Nkx,1:Nz), ionode=0) ELSE ! output using putarrnd (very slow on marconi) CALL putarrnd(fidres, dset_name, field(ikys:ikye,ikxs:ikxe, izs:ize), (/1, 1, 3/)) ENDIF CALL attach(fidres, dset_name, "time", time) END SUBROUTINE write_field3d_kykxz SUBROUTINE write_field3d_pjz_i(field, text) IMPLICIT NONE REAL(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize), INTENT(IN) :: field CHARACTER(*), INTENT(IN) :: text CHARACTER(LEN=50) :: dset_name WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d IF (num_procs .EQ. 1) THEN ! no data distribution CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,izs:ize), ionode=0) ELSE CALL putarrnd(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,izs:ize), (/1, 0, 3/)) ENDIF CALL attach(fidres, dset_name, "time", time) END SUBROUTINE write_field3d_pjz_i SUBROUTINE write_field3d_pjz_e(field, text) IMPLICIT NONE REAL(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize), INTENT(IN) :: field CHARACTER(*), INTENT(IN) :: text CHARACTER(LEN=50) :: dset_name WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d IF (num_procs .EQ. 1) THEN ! no data distribution CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,izs:ize), ionode=0) ELSE CALL putarrnd(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,izs:ize), (/1, 0, 3/)) ENDIF CALL attach(fidres, dset_name, "time", time) END SUBROUTINE write_field3d_pjz_e END SUBROUTINE diagnose_3d SUBROUTINE diagnose_5d USE basic USE futils, ONLY: append, getatt, attach, putarrnd, putarr 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, & Np_i, Nj_i, Np_e, Nj_e, Nky, Nkx, Nz, & ikxs,ikxe,ikys,ikye,izs,ize USE time_integration USE diagnostics_par USE prec_const USE model, ONLY: KIN_E IMPLICIT NONE CHARACTER(LEN=50) :: dset_name 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) IF (write_Napj) THEN IF(KIN_E)& CALL write_field5d_e(moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel), 'moments_e') CALL write_field5d_i(moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel), 'moments_i') ENDIF IF (write_Sapj) THEN IF(KIN_E)& CALL write_field5d_e(Sepj(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), 'Sepj') CALL write_field5d_i(Sipj(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), 'Sipj') ENDIF CONTAINS SUBROUTINE write_field5d_e(field, text) USE futils, ONLY: attach, putarr, putarrnd USE parallel, ONLY: gather_pjxyz_e USE grid, ONLY: ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize USE prec_const IMPLICIT NONE COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), INTENT(IN) :: field CHARACTER(*), INTENT(IN) :: text COMPLEX(dp), DIMENSION(1:Np_e,1:Nj_e,1:Nky,1:Nkx,1:Nz) :: field_full 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,ikys:ikye,ikxs:ikxe,izs:ize), ionode=0) ELSEIF(GATHERV_OUTPUT) THEN ! output using one node (gatherv) CALL gather_pjxyz_e(field(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize),& field_full(1:Np_e,1:Nj_e,1:Nky,1:Nkx,1:Nz)) CALL putarr(fidres, dset_name, field_full(1:Np_i,1:Nj_i,1:Nky,1:Nkx,1:Nz), ionode=0) ELSE CALL putarrnd(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), (/1,3,5/)) 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 parallel, ONLY: gather_pjxyz_i USE grid, ONLY: ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize USE prec_const IMPLICIT NONE COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), INTENT(IN) :: field CHARACTER(*), INTENT(IN) :: text COMPLEX(dp), DIMENSION(1:Np_i,1:Nj_i,1:Nky,1:Nkx,1:Nz) :: field_full 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,ikys:ikye,ikxs:ikxe,izs:ize), ionode=0) ELSEIF(GATHERV_OUTPUT) THEN ! output using one node (gatherv) CALL gather_pjxyz_i(field(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize),& field_full(1:Np_i,1:Nj_i,1:Nky,1:Nkx,1:Nz)) CALL putarr(fidres, dset_name, field_full(1:Np_i,1:Nj_i,1:Nky,1:Nkx,1:Nz), ionode=0) ELSE CALL putarrnd(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), (/1,3,5/)) 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 SUBROUTINE spit_snapshot_check USE fields, ONLY: phi USE grid, ONLY: ikxs,ikxe,Nkx,ikys,ikye,Nky,izs,ize,Nz USE parallel, ONLY: gather_xyz USE basic IMPLICIT NONE LOGICAL :: file_exist INTEGER :: fid_check, ikx, iky, iz CHARACTER(256) :: check_filename COMPLEX(dp), DIMENSION(1:Nky,1:Nkx,1:Nz) :: field_to_check !! Spit a snapshot of PHI if requested (triggered by creating a file named "check_phi") INQUIRE(file='check_phi', exist=file_exist) IF( file_exist ) THEN IF(my_id.EQ. 0) WRITE(*,*) 'Check file found -> gather phi..' CALL gather_xyz(phi(ikys:ikye,ikxs:ikxe,izs:ize), field_to_check) IF(my_id.EQ. 0) THEN WRITE(check_filename,'(a16)') 'check_phi.out' OPEN(fid_check, file=check_filename, form='formatted') WRITE(*,*) 'Check file found -> output phi ..' WRITE(fid_check,*) Nky, Nkx, Nz DO iky = 1,Nky; DO ikx = 1, Nkx; DO iz = 1,Nz WRITE(fid_check,*) real(field_to_check(iky,ikx,iz)), ',' , imag(field_to_check(iky,ikx,iz)) ENDDO; ENDDO; ENDDO CLOSE(fid_check) WRITE(*,*) 'Check file found -> done.' ! delete the check_phi flagfile OPEN(fid_check, file='check_phi') CLOSE(fid_check, status='delete') ENDIF ENDIF END SUBROUTINE spit_snapshot_check diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90 index 9840cc0..4921117 100644 --- a/src/geometry_mod.F90 +++ b/src/geometry_mod.F90 @@ -1,592 +1,432 @@ module geometry ! computes geometrical quantities ! Adapted from B.J.Frei MOLIX code (2021) use prec_const use model use grid use array use fields use basic use calculus, ONLY: simpson_rule_z use miller, ONLY: set_miller_parameters, get_miller implicit none PRIVATE ! Geometry input parameters CHARACTER(len=16), & PUBLIC, PROTECTED :: geom REAL(dp), PUBLIC, PROTECTED :: q0 = 1.4_dp ! safety factor REAL(dp), PUBLIC, PROTECTED :: shear = 0._dp ! magnetic field shear REAL(dp), PUBLIC, PROTECTED :: eps = 0.18_dp ! inverse aspect ratio REAL(dp), PUBLIC, PROTECTED :: alpha_MHD = 0 ! shafranov shift effect alpha = -q2 R dbeta/dr ! parameters for Miller geometry REAL(dp), PUBLIC, PROTECTED :: kappa = 1._dp ! elongation REAL(dp), PUBLIC, PROTECTED :: s_kappa = 0._dp ! r normalized derivative skappa = r/kappa dkappa/dr REAL(dp), PUBLIC, PROTECTED :: delta = 0._dp ! triangularity REAL(dp), PUBLIC, PROTECTED :: s_delta = 0._dp ! '' sdelta = r/sqrt(1-delta2) ddelta/dr REAL(dp), PUBLIC, PROTECTED :: zeta = 0._dp ! squareness REAL(dp), PUBLIC, PROTECTED :: s_zeta = 0._dp ! '' szeta = r dzeta/dr ! GENE unused additional parameters for miller_mod REAL(dp), PUBLIC, PROTECTED :: edge_opt = 0 ! meant to redistribute the points in z REAL(dp), PUBLIC, PROTECTED :: major_R = 1 ! major radius REAL(dp), PUBLIC, PROTECTED :: major_Z = 0 ! vertical elevation REAL(dp), PUBLIC, PROTECTED :: dpdx_pm_geom = 0 ! amplitude mag. eq. pressure grad. REAL(dp), PUBLIC, PROTECTED :: C_y = 0 ! defines y coordinate : Cy (q theta - phi) REAL(dp), PUBLIC, PROTECTED :: C_xy = 0 ! defines x coordinate : B = Cxy Vx x Vy ! Geometrical auxiliary variables LOGICAL, PUBLIC, PROTECTED :: SHEARED = .false. ! flag for shear magn. geom or not ! Curvature REAL(dp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z, odd/even p ! Jacobian REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Jacobian ! dimensions: z, odd/even p COMPLEX(dp), PUBLIC, PROTECTED :: iInt_Jacobian ! Inverse integrated Jacobian ! Metric REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gxx, gxy, gxz, gyy, gyz, gzz ! dimensions: z, odd/even p REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dxdr, dxdZ, Rc, phic, Zc ! derivatives of magnetic field strength REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gradxB, gradyB, gradzB ! Relative magnetic field strength REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB, hatB_NL ! Relative strength of major radius REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ ! Some geometrical coefficients REAL(dp), PUBLIC, DIMENSION(:,:) , ALLOCATABLE :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ] ! Array to map the index of mode (kx,ky,-pi) to (kx+2pi*s*ky,ky,pi) for sheared periodic boundary condition INTEGER, PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ikx_zBC_L, ikx_zBC_R ! Functions PUBLIC :: geometry_readinputs, geometry_outputinputs,& eval_magnetic_geometry, set_ikx_zBC_map CONTAINS SUBROUTINE geometry_readinputs ! Read the input parameters IMPLICIT NONE NAMELIST /GEOMETRY/ geom, q0, shear, eps,& kappa, s_kappa,delta, s_delta, zeta, s_zeta ! For miller READ(lu_in,geometry) IF(shear .NE. 0._dp) SHEARED = .true. END SUBROUTINE geometry_readinputs subroutine eval_magnetic_geometry ! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient implicit none REAL(dp) :: kx,ky COMPLEX(dp), DIMENSION(izs:ize) :: integrant INTEGER :: fid ! Allocate arrays CALL geometry_allocate_mem ! IF( (Ny .EQ. 1) .AND. (Nz .EQ. 1)) THEN !1D perp linear run IF( my_id .eq. 0 ) WRITE(*,*) '1D perpendicular geometry' call eval_1D_geometry ELSE SELECT CASE(geom) - CASE('circular') - IF( my_id .eq. 0 ) WRITE(*,*) 'circular geometry' - call eval_circular_geometry CASE('s-alpha') IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry' call eval_salphaB_geometry CASE('Z-pinch') IF( my_id .eq. 0 ) WRITE(*,*) 'Z-pinch geometry' call eval_zpinch_geometry SHEARED = .FALSE. CASE('miller') IF( my_id .eq. 0 ) WRITE(*,*) 'Miller geometry' call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta) call get_miller(eps,major_R,major_Z,q0,shear,alpha_MHD,edge_opt,& C_y,C_xy,dpdx_pm_geom,gxx,gyy,gzz,gxy,gxz,gyz,& - gradxB,gradyB,hatB,jacobian,gradzB,hatR,hatZ,dxdR,dxdZ) + gradxB,gradyB,hatB,jacobian,gradzB,hatR,hatZ,dxdR,dxdZ,& + Ckxky,gradz_coeff) CASE DEFAULT STOP 'geometry not recognized!!' END SELECT ENDIF ! ! Evaluate perpendicular wavenumber ! k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy} ! normalized to rhos_ DO eo = 0,1 - DO iky = ikys, ikye - ky = kyarray(iky) - DO ikx = ikxs, ikxe - kx = kxarray(ikx) - DO iz = izgs,izge - kparray(iky, ikx, iz, eo) = & - SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/hatB(iz,eo) - ! there is a factor 1/B from the normalization; important to match GENE + DO iz = izgs,izge + DO iky = ikys, ikye + ky = kyarray(iky) + DO ikx = ikxs, ikxe + kx = kxarray(ikx) + kparray(iky, ikx, iz, eo) = & + SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/hatB(iz,eo) + ! there is a factor 1/B from the normalization; important to match GENE ENDDO - ENDDO - ENDDO + ENDDO + ENDDO ENDDO + ! Factor in front of the nonlinear term + hatB_NL(izgs:izge,0:1) = Jacobian(izgs:izge,0:1)& + *(gxx(izgs:izge,0:1)*gyy(izgs:izge,0:1) - gxy(izgs:izge,0:1)**2)/hatB(izgs:izge,0:1) + ! set the mapping for parallel boundary conditions CALL set_ikx_zBC_map two_third_kpmax = 2._dp/3._dp * MAXVAL(kparray) ! ! Compute the inverse z integrated Jacobian (useful for flux averaging) integrant = Jacobian(izs:ize,0) ! Convert into complex array CALL simpson_rule_z(integrant,iInt_Jacobian) iInt_Jacobian = 1._dp/iInt_Jacobian ! reverse it END SUBROUTINE eval_magnetic_geometry ! !-------------------------------------------------------------------------------- ! SUBROUTINE eval_salphaB_geometry ! evaluate s-alpha geometry model implicit none REAL(dp) :: z, kx, ky, Gx, Gy alpha_MHD = 0._dp parity: DO eo = 0,1 zloop: DO iz = izgs,izge z = zarray(iz,eo) ! metric gxx(iz,eo) = 1._dp gxy(iz,eo) = shear*z - alpha_MHD*SIN(z) gxz(iz,eo) = 0._dp gyy(iz,eo) = 1._dp + (shear*z - alpha_MHD*SIN(z))**2 gyz(iz,eo) = 1._dp/eps gzz(iz,eo) = 0._dp dxdR(iz,eo)= COS(z) dxdZ(iz,eo)= SIN(z) ! Relative strengh of radius hatR(iz,eo) = 1._dp + eps*COS(z) hatZ(iz,eo) = 1._dp + eps*SIN(z) ! toroidal coordinates Rc (iz,eo) = hatR(iz,eo) phic(iz,eo) = z Zc (iz,eo) = hatZ(iz,eo) ! Jacobian Jacobian(iz,eo) = q0*hatR(iz,eo) ! Relative strengh of modulus of B - hatB (iz,eo) = 1._dp / hatR(iz,eo) - hatB_NL(iz,eo) = 1._dp ! Factor in front of the nonlinear term + hatB(iz,eo) = 1._dp / hatR(iz,eo) ! Derivative of the magnetic field strenght gradxB(iz,eo) = -COS(z) ! Gene put a factor hatB^2 or 1/hatR^2 in this gradyB(iz,eo) = 0._dp gradzB(iz,eo) = eps*SIN(z)/hatR(iz,eo) ! Gene put a factor hatB or 1/hatR in this ! Curvature operator Gx = (gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo)) *eps*SIN(Z) ! Kx Gy = -COS(z) + (gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)) *eps*SIN(Z) ! Ky DO iky = ikys, ikye ky = kyarray(iky) DO ikx= ikxs, ikxe kx = kxarray(ikx) - Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - COS(z)*ky -(shear*z-alpha_MHD*SIN(z))*SIN(z)*ky)* hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine + Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - COS(z)*ky -(shear*z-alpha_MHD*SIN(z))*SIN(z)*ky)/ hatB(iz,eo) ! Ckxky(iky, ikx, iz,eo) = (Gx*kx + Gy*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine ENDDO ENDDO ! coefficient in the front of parallel derivative gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo) ENDDO zloop ENDDO parity END SUBROUTINE eval_salphaB_geometry ! !-------------------------------------------------------------------------------- ! - - SUBROUTINE eval_circular_geometry - ! evaluate circular geometry model - ! Ref: Lapilonne et al., PoP, 2009, GENE circular.F90 applied at r=r0 + SUBROUTINE eval_zpinch_geometry implicit none - REAL(dp) :: chi, kx, ky, Gx, Gy + REAL(dp) :: z, kx, ky, alpha_MHD + alpha_MHD = 0._dp - parity: DO eo = 0,1 - zloop: DO iz = izgs,izge - chi = zarray(iz,eo) ! = chi + parity: DO eo = 0,1 + zloop: DO iz = izgs,izge + z = zarray(iz,eo) - ! metric in x,y,z + ! metric gxx(iz,eo) = 1._dp - gyy(iz,eo) = 1._dp + (shear*chi)**2 - 2._dp*eps*COS(chi) - 2._dp*shear*chi*eps*SIN(chi) - gxy(iz,eo) = shear*chi - eps*SIN(chi) - gxz(iz,eo) = -SIN(chi) - gyz(iz,eo) = (1._dp - 2._dp*eps*COS(chi) - shear*chi*eps*SIN(chi))/eps - gzz(iz,eo) = (1._dp - 2._dp*eps*COS(chi))/eps**2 - dxdR(iz,eo)= COS(chi) - dxdZ(iz,eo)= SIN(chi) + gxy(iz,eo) = 0._dp + gxz(iz,eo) = 0._dp + gyy(iz,eo) = 1._dp + gyz(iz,eo) = 0._dp + gzz(iz,eo) = 1._dp + dxdR(iz,eo)= COS(z) + dxdZ(iz,eo)= SIN(z) ! Relative strengh of radius - hatR(iz,eo) = 1._dp + eps*COS(chi) - hatZ(iz,eo) = 1._dp + eps*SIN(chi) + hatR(iz,eo) = 1._dp + hatZ(iz,eo) = 1._dp ! toroidal coordinates Rc (iz,eo) = hatR(iz,eo) - phic(iz,eo) = chi + phic(iz,eo) = z Zc (iz,eo) = hatZ(iz,eo) ! Jacobian - Jacobian(iz,eo) = q0*hatR(iz,eo) + Jacobian(iz,eo) = 1._dp ! Relative strengh of modulus of B - hatB (iz,eo) = 1._dp / hatR(iz,eo) - hatB_NL(iz,eo) = 1._dp ! Factor in front of the nonlinear term + hatB (iz,eo) = 1._dp + hatB_NL(iz,eo) = 1._dp ! Derivative of the magnetic field strenght - gradxB(iz,eo) = -COS(chi) ! Gene put a factor hatB^2 or 1/hatR^2 in this - gradyB(iz,eo) = 0._dp - gradzB(iz,eo) = eps * SIN(chi) / hatR(iz,eo) ! Gene put a factor hatB or 1/hatR in this + gradxB(iz,eo) = 0._dp ! Gene put a factor hatB^2 or 1/hatR^2 in this + gradyB(iz,eo) = 0._dp + gradzB(iz,eo) = 0._dp ! Gene put a factor hatB or 1/hatR in this ! Curvature operator - Gx = (gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo)) *eps*SIN(chi) ! Kx - Gy = -COS(chi) + (gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)) *eps*SIN(chi) ! Ky - DO iky = ikys, ikye + DO iky = ikys, ikye ky = kyarray(iky) DO ikx= ikxs, ikxe kx = kxarray(ikx) - Ckxky(iky, ikx, iz,eo) = (Gx*kx + Gy*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine + Ckxky(iky, ikx, iz,eo) = -ky * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine ENDDO ENDDO ! coefficient in the front of parallel derivative gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo) - ENDDO zloop - ENDDO parity - - END SUBROUTINE eval_circular_geometry - ! - !-------------------------------------------------------------------------------- - ! - - SUBROUTINE eval_circular_geometry_GENE - ! evaluate circular geometry model - ! Ref: Lapilonne et al., PoP, 2009, GENE circular.F90 applied at r=r0 - implicit none - REAL(dp) :: qbar, dxdr_, dpsidr, dqdr, dqbardr, Cy, dCydr_Cy, Cxy, fac2, & - cost, sint, dchidr, dchidt, sign_Ip, B, dBdt, dBdr, & - g11, g22, g12, g33, dBdchi, dBdr_c !GENE variables - REAL(dp) :: X, z, kx, ky, Gamma1, Gamma2, Gamma3, feps - - sign_Ip = 1._dp - feps = SQRT(1._dp-eps**2) - qbar = q0 * feps - dxdr_ = 1._dp - dpsidr = eps/qbar - dqdr = q0/eps*shear - dqbardr = feps*q0/eps*(shear - eps**2/feps**2) - Cy = eps/q0 * sign_Ip - dCydr_Cy = 0._dp - Cxy = 1._dp/feps - fac2 = dCydr_Cy + dqdr/q0 - - - parity: DO eo = 0,1 - zloop: DO iz = izgs,izge - z = zarray(iz,eo) - - cost = (COS(z)-eps)/(1._dp-eps*COS(z)) - sint = feps*sin(sign_Ip*z)/(1._dp-eps*COS(z)) - dchidr = -SIN(z)/feps**2 - dchidt = sign_Ip*feps/(1._dp+eps*cost) - B = SQRT(1._dp+(eps/qbar)**2)/(1._dp+eps*cost) - dBdt = eps*sint*B/(1._dp+eps*cost) - - ! metric in r,chi,Phi - g11 = 1._dp - g22 = dchidr**2 + dchidt**2/eps**2 - g12 = dchidr - g33 = 1._dp/(1._dp+eps*cost)**2 - - ! magnetic field derivatives in - dBdchi=dBdt/dchidt - dBdr_c=(dBdr-g12*dBdt/dchidt)/g11 - - ! metric in x,y,z - gxx(iz,eo) = dxdr_**2*g11 - gyy(iz,eo) = (Cy*q0)**2 * (fac2*z)**2*g11 & - + 2._dp*fac2*z*g12 & - + Cy**2*g33 - gxy(iz,eo) = dxdr_*Cy*sign_Ip*q0*(fac2*z*g11 + g12) - gxz(iz,eo) = dxdr_*g12 - gyz(iz,eo) = Cy*q0*sign_Ip*(fac2*z*g12 + g22) - gzz(iz,eo) = g22 - - ! Jacobian - Jacobian(iz,eo) = Cxy*abs(q0)*(1._dp+eps*cost)**2 - - ! Background equilibrium magnetic field - hatB (iz,eo) = B - hatB_NL(iz,eo) = SQRT(gxx(iz,eo)*gyy(iz,eo) - (gxy(iz,eo))**2) ! In front of the NL term - gradxB(iz,eo) = dBdr_c/dxdr_ ! Gene put a factor hatB^2 or 1/hatR^2 in this - gradyB(iz,eo) = 0._dp - gradzB(iz,eo) = dBdchi! Gene put a factor hatB or 1/hatR in this - - ! Cylindrical coordinates derivatives - dxdR(iz,eo) = cost - dxdZ(iz,eo) = sint - hatR(iz,eo) = 1._dp + eps*cost - hatZ(iz,eo) = 1._dp + eps*sint - - ! toroidal coordinates - Rc (iz,eo) = hatR(iz,eo) - phic(iz,eo) = X - Zc (iz,eo) = hatZ(iz,eo) - - Gamma1 = gxy(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyy(iz,eo) - Gamma2 = gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo) ! Kx - Gamma3 = gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo) ! Ky - ! Curvature operator - DO iky = ikys, ikye - ky = kyarray(iky) - DO ikx= ikxs, ikxe - kx = kxarray(ikx) - ! Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - (COS(z) + (shear*z - alpha_MHD*SIN(z))* SIN(z))*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine - Ckxky(iky, ikx, iz,eo) = (-sint*kx - cost*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine - ENDDO - ENDDO - ! coefficient in the front of parallel derivative - gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo) - - ENDDO zloop - ENDDO parity - - END SUBROUTINE eval_circular_geometry_GENE - ! - !-------------------------------------------------------------------------------- - ! - - SUBROUTINE eval_zpinch_geometry - ! evaluate s-alpha geometry model - implicit none - REAL(dp) :: z, kx, ky, alpha_MHD - alpha_MHD = 0._dp - - parity: DO eo = 0,1 - zloop: DO iz = izgs,izge - z = zarray(iz,eo) - - ! metric - gxx(iz,eo) = 1._dp - gxy(iz,eo) = 0._dp - gxz(iz,eo) = 0._dp - gyy(iz,eo) = 1._dp - gyz(iz,eo) = 0._dp - gzz(iz,eo) = 1._dp - dxdR(iz,eo)= COS(z) - dxdZ(iz,eo)= SIN(z) - - ! Relative strengh of radius - hatR(iz,eo) = 1._dp - hatZ(iz,eo) = 1._dp - - ! toroidal coordinates - Rc (iz,eo) = hatR(iz,eo) - phic(iz,eo) = z - Zc (iz,eo) = hatZ(iz,eo) - - ! Jacobian - Jacobian(iz,eo) = 1._dp - - ! Relative strengh of modulus of B - hatB (iz,eo) = 1._dp - hatB_NL(iz,eo) = 1._dp - - ! Derivative of the magnetic field strenght - gradxB(iz,eo) = 0._dp ! Gene put a factor hatB^2 or 1/hatR^2 in this - gradyB(iz,eo) = 0._dp - gradzB(iz,eo) = 0._dp ! Gene put a factor hatB or 1/hatR in this - - ! Curvature operator - DO iky = ikys, ikye - ky = kyarray(iky) - DO ikx= ikxs, ikxe - kx = kxarray(ikx) - Ckxky(iky, ikx, iz,eo) = -ky * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine - ENDDO - ENDDO - ! coefficient in the front of parallel derivative - gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo) - - ENDDO zloop - ENDDO parity + ENDDO zloop + ENDDO parity END SUBROUTINE eval_zpinch_geometry ! !-------------------------------------------------------------------------------- ! subroutine eval_1D_geometry ! evaluate 1D perp geometry model implicit none REAL(dp) :: z, kx, ky parity: DO eo = 0,1 zloop: DO iz = izs,ize z = zarray(iz,eo) ! metric gxx(iz,eo) = 1._dp gxy(iz,eo) = 0._dp gyy(iz,eo) = 1._dp ! Relative strengh of radius hatR(iz,eo) = 1._dp ! Jacobian Jacobian(iz,eo) = 1._dp ! Relative strengh of modulus of B hatB(iz,eo) = 1._dp ! Curvature operator DO iky = ikys, ikye ky = kyarray(iky) DO ikx= ikxs, ikxe kx = kxarray(ikx) Ckxky(ikx, iky, iz,eo) = -kx ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine ENDDO ENDDO ! coefficient in the front of parallel derivative gradz_coeff(iz,eo) = 1._dp ENDDO zloop ENDDO parity END SUBROUTINE eval_1D_geometry ! !-------------------------------------------------------------------------------- ! SUBROUTINE set_ikx_zBC_map IMPLICIT NONE REAL :: shift, kx_shift ! For periodic CHI BC or 0 dirichlet LOGICAL :: PERIODIC_CHI_BC = .TRUE. ALLOCATE(ikx_zBC_R(ikys:ikye,ikxs:ikxe)) ALLOCATE(ikx_zBC_L(ikys:ikye,ikxs:ikxe)) ! No periodic connection for extension of the domain IF(Nexc .GT. 1) PERIODIC_CHI_BC = .TRUE. !! No shear case (simple id mapping) !3 | 1 2 3 4 5 6 | ky = 3 dky !2 ky | 1 2 3 4 5 6 | ky = 2 dky !1 A | 1 2 3 4 5 6 | ky = 1 dky !0 | -> kx | 1____2____3____4____5____6 | ky = 0 dky !(e.g.) kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=free) DO iky = ikys,ikye DO ikx = ikxs,ikxe ikx_zBC_L(iky,ikx) = ikx ikx_zBC_R(iky,ikx) = ikx ENDDO ENDDO ! Modify connection map only at border of z IF(SHEARED) THEN ! connection map BC of the RIGHT boundary (z=pi*Npol-dz) (even NZ) !3 | 4 x x x 2 3 | ky = 3 dky !2 ky | 3 4 x x 1 2 | ky = 2 dky !1 A | 2 3 4 x 6 1 | ky = 1 dky !0 | -> kx | 1____2____3____4____5____6 | ky = 0 dky !kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=2pi*shear*npol*dky) ! connection map BC of the RIGHT boundary (z=pi*Npol-dz) (ODD NZ) !3 | x x x 2 3 | ky = 3 dky !2 ky | 3 x x 1 2 | ky = 2 dky !1 A | 2 3 x 5 1 | ky = 1 dky !0 | -> kx | 1____2____3____4____5 | ky = 0 dky !kx = 0 0.1 0.2 -0.2 -0.1 (dkx=2pi*shear*npol*dky) IF(contains_zmax) THEN ! Check if the process is at the end of the FT DO iky = ikys,ikye shift = 2._dp*PI*shear*kyarray(iky)*Npol DO ikx = ikxs,ikxe kx_shift = kxarray(ikx) + shift ! We use EPSILON() to treat perfect equality case IF( ((kx_shift-EPSILON(kx_shift)) .GT. kx_max) .AND. (.NOT. PERIODIC_CHI_BC) )THEN ! outside of the frequ domain ikx_zBC_R(iky,ikx) = -99 ELSE ikx_zBC_R(iky,ikx) = ikx+(iky-1)*Nexc IF( ikx_zBC_R(iky,ikx) .GT. Nkx ) & ikx_zBC_R(iky,ikx) = ikx_zBC_R(iky,ikx) - Nkx ENDIF ENDDO ENDDO ENDIF ! connection map BC of the LEFT boundary (z=-pi*Npol) !3 | x 5 6 1 x x | ky = 3 dky !2 ky | 5 6 1 2 x x | ky = 2 dky !1 A | 6 1 2 3 x 5 | ky = 1 dky !0 | -> kx | 1____2____3____4____5____6 | ky = 0 dky !(e.g.) kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=2pi*shear*npol*dky) IF(contains_zmin) THEN ! Check if the process is at the start of the FT DO iky = ikys,ikye shift = 2._dp*PI*shear*kyarray(iky)*Npol DO ikx = ikxs,ikxe kx_shift = kxarray(ikx) - shift ! We use EPSILON() to treat perfect equality case IF( ((kx_shift+EPSILON(kx_shift)) .LT. kx_min) .AND. (.NOT. PERIODIC_CHI_BC) ) THEN ! outside of the frequ domain ikx_zBC_L(iky,ikx) = -99 ELSE ikx_zBC_L(iky,ikx) = ikx-(iky-1)*Nexc IF( ikx_zBC_L(iky,ikx) .LT. 1 ) & ikx_zBC_L(iky,ikx) = ikx_zBC_L(iky,ikx) + Nkx ENDIF ENDDO ENDDO ENDIF ELSE ENDIF END SUBROUTINE set_ikx_zBC_map ! !-------------------------------------------------------------------------------- ! SUBROUTINE geometry_allocate_mem ! Curvature and geometry CALL allocate_array( Ckxky, ikys,ikye, ikxs,ikxe,izgs,izge,0,1) CALL allocate_array( Jacobian,izgs,izge, 0,1) CALL allocate_array( gxx,izgs,izge, 0,1) CALL allocate_array( gxy,izgs,izge, 0,1) CALL allocate_array( gxz,izgs,izge, 0,1) CALL allocate_array( gyy,izgs,izge, 0,1) CALL allocate_array( gyz,izgs,izge, 0,1) CALL allocate_array( gzz,izgs,izge, 0,1) CALL allocate_array( gradxB,izgs,izge, 0,1) CALL allocate_array( gradyB,izgs,izge, 0,1) CALL allocate_array( gradzB,izgs,izge, 0,1) CALL allocate_array( hatB,izgs,izge, 0,1) CALL allocate_array( hatB_NL,izgs,izge, 0,1) CALL allocate_array( hatR,izgs,izge, 0,1) CALL allocate_array( hatZ,izgs,izge, 0,1) CALL allocate_array( Rc,izgs,izge, 0,1) CALL allocate_array( phic,izgs,izge, 0,1) CALL allocate_array( Zc,izgs,izge, 0,1) CALL allocate_array( dxdR,izgs,izge, 0,1) CALL allocate_array( dxdZ,izgs,izge, 0,1) call allocate_array(gradz_coeff,izgs,izge, 0,1) CALL allocate_array( kparray, ikys,ikye, ikxs,ikxe,izgs,izge,0,1) END SUBROUTINE geometry_allocate_mem SUBROUTINE geometry_outputinputs(fidres, str) ! Write the input parameters to the results_xx.h5 file USE futils, ONLY: attach USE prec_const IMPLICIT NONE INTEGER, INTENT(in) :: fidres CHARACTER(len=256), INTENT(in) :: str CALL attach(fidres, TRIM(str),"geometry", geom) CALL attach(fidres, TRIM(str), "q0", q0) CALL attach(fidres, TRIM(str), "shear", shear) CALL attach(fidres, TRIM(str), "eps", eps) END SUBROUTINE geometry_outputinputs end module geometry diff --git a/src/miller_mod.F90 b/src/miller_mod.F90 index f75a075..b33e3a5 100644 --- a/src/miller_mod.F90 +++ b/src/miller_mod.F90 @@ -1,645 +1,667 @@ !! This source has been adapted from GENE https://genecode.org/ !! !>Implementation of the local equilibrium model of [R.L. Miller et al., PoP 5, 973 (1998) !>and [J. Candy, PPCF 51, 105009 (2009)] MODULE miller USE prec_const USE basic ! use coordinates,only: gcoor, get_dzprimedz USE grid ! use discretization USE lagrange_interpolation ! use par_in, only: beta, sign_Ip_CW, sign_Bt_CW, Npol USE model implicit none public:: get_miller, set_miller_parameters public:: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta public:: thetaShift public:: mMode, nMode public:: thetak, thetad public:: aSurf, Delta2, Delta3, theta2, theta3, Raxis, Zaxis public:: Deltam, Deltan, s_Deltam, s_Deltan, thetam, thetan public:: cN_m, sN_m, cNdr_m, sNdr_m private real(dp) :: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta INTEGER :: mMode, nMode real(dp) :: thetaShift real(dp) :: thetak, thetad real(dp) :: aSurf, Delta2, Delta3, theta2, theta3, Raxis, Zaxis real(dp) :: Deltam, Deltan, s_Deltam, s_Deltan, thetam, thetan INTEGER, PARAMETER :: IND_M=32 real(dp), DIMENSION(0:IND_M-1) :: cN_m, sN_m, cNdr_m, sNdr_m CONTAINS !>Set defaults for miller parameters subroutine set_miller_parameters(kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_) real(dp), INTENT(IN) :: kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_ rho = -1.0 kappa = kappa_ s_kappa = s_kappa_ delta = delta_ s_delta = s_delta_ zeta = zeta_ s_zeta = s_zeta_ drR = 0.0 drZ = 0.0 thetak = 0.0 thetad = 0.0 aSurf = 0.54 Delta2 = 1.0 Delta3 = 1.0 theta2 = 0.0 theta3 = 0.0 Raxis = 1.0 Zaxis = 0.0 mMode = 2 nMode = 3 Deltam = 1.0 Deltan = 1.0 s_Deltam = 0.0 s_Deltan = 0.0 thetam = 0.0 thetan = 0.0 cN_m = 0.0 sN_m = 0.0 cNdr_m = 0.0 sNdr_m = 0.0 end subroutine set_miller_parameters !>Get Miller metric, magnetic field, jacobian etc. subroutine get_miller(trpeps,major_R,major_Z,q0,shat,amhd,edge_opt,& C_y,C_xy,dpdx_pm_geom,gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,dBdx_,dBdy_,& - Bfield_,jacobian_,dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_) + Bfield_,jacobian_,dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_,Ckxky_,gradz_coeff_) !!!!!!!!!!!!!!!! GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real(dp), INTENT(INOUT) :: trpeps ! eps in gyacomo (inverse aspect ratio) real(dp), INTENT(INOUT) :: major_R ! major radius real(dp), INTENT(INOUT) :: major_Z ! major Z real(dp), INTENT(INOUT) :: q0 ! safetyfactor real(dp), INTENT(INOUT) :: shat ! safetyfactor real(dp), INTENT(INOUT) :: amhd ! alpha mhd real(dp), INTENT(INOUT) :: edge_opt ! alpha mhd real(dp), INTENT(INOUT) :: dpdx_pm_geom ! amplitude mag. eq. pressure grad. real(dp), INTENT(INOUT) :: C_y, C_xy real(dp), dimension(izgs:izge,0:1), INTENT(INOUT) :: & gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,& dBdx_,dBdy_,Bfield_,jacobian_,& - dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_ + dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_, & + gradz_coeff_ + real(dp), dimension(ikys:ikye,ikxs:ikxe,izgs:izge,0:1), INTENT(INOUT) :: Ckxky_ ! No parameter in gyacomo yet integer :: ikxgs =1 ! the left ghost gpdisc%pi1gl real(dp) :: sign_Ip_CW=1 ! current sign (only normal current) real(dp) :: sign_Bt_CW=1 ! current sign (only normal current) !!!!!!!!!!!!!! END GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - + ! Auxiliary variables for curvature computation + real(dp) :: G1,G2,G3,Cx,Cy,ky,kx integer:: np, np_s, Npol_ext, Npol_s real(dp), dimension(500*(Npol+2)):: R,Z,R_rho,Z_rho,R_theta,Z_theta,R_theta_theta,Z_theta_theta,dlp,Rc,cosu,sinu,Bphi real(dp), dimension(500*(Npol+2)):: drRcirc, drRelong, drRelongTilt, drRtri, drRtriTilt, drZcirc, drZelong, drZelongTilt real(dp), dimension(500*(Npol+2)):: drZtri, drZtriTilt, dtdtRcirc, dtdtRelong, dtdtRelongTilt, dtdtRtri, dtdtRtriTilt real(dp), dimension(500*(Npol+2)):: dtdtZcirc, dtdtZelong, dtdtZelongTilt, dtdtZtri, dtdtZtriTilt, dtRcirc, dtRelong real(dp), dimension(500*(Npol+2)):: dtRelongTilt, dtRtri, dtRtriTilt, dtZcirc, dtZelong, dtZelongTilt, dtZtri, dtZtriTilt real(dp), dimension(500*(Npol+2)):: Rcirc, Relong, RelongTilt, Rtri, RtriTilt, Zcirc, Zelong, ZelongTilt, Ztri, ZtriTilt real(dp), dimension(500*(Npol+2)):: drrShape, drrAng, drxAng, dryAng, dtdtrShape, dtdtrAng, dtdtxAng real(dp), dimension(500*(Npol+2)):: dtdtyAng, dtrShape, dtrAng, dtxAng, dtyAng, rShape, rAng, xAng, yAng real(dp), dimension(500*(Npol+2)):: theta, thAdj, J_r, B, Bp, D0, D1, D2, D3, nu, chi, psi1, nu1 real(dp), dimension(500*(Npol+2)):: tmp_reverse, theta_reverse, tmp_arr real(dp), dimension(500*(Npol+1)):: theta_s, thAdj_s, chi_s, theta_s_reverse real(dp), dimension(500*(Npol+1)):: R_s, Z_s, R_theta_s, Z_theta_s, Rc_s, cosu_s, sinu_s, Bphi_s, B_s, Bp_s, dlp_s real(dp), dimension(500*(Npol+1)):: dtRcirc_s, dtRelong_s, dtRelongTilt_s, dtRtri_s, dtRtriTilt_s, dtZcirc_s real(dp), dimension(500*(Npol+1)):: dtZelong_s, dtZelongTilt_s, dtZtri_s, dtZtriTilt_s, Rcirc_s, Relong_s, RelongTilt_s real(dp), dimension(500*(Npol+1)):: Rtri_s, RtriTilt_s, Zcirc_s, Zelong_s, ZelongTilt_s, Ztri_s, ZtriTilt_s, dtrShape_s real(dp), dimension(500*(Npol+1)):: dtrAng_s, dtxAng_s, dtyAng_s, rShape_s, rAng_s, xAng_s, yAng_s real(dp), dimension(500*(Npol+1)):: psi1_s, nu1_s, dchidx_s, dB_drho_s, dB_dl_s, dnu_drho_s, dnu_dl_s, grad_nu_s real(dp), dimension(500*(Npol+1)):: gxx, gxy, gxz, gyy, gyz, gzz, dtheta_dchi_s, dBp_dchi_s, jacobian, dBdx, dBdz real(dp), dimension(500*(Npol+1)):: g_xx, g_xy, g_xz, g_yy, g_yz, g_zz, tmp_arr_s, dxdR_s, dxdZ_s, K_x, K_y !tmp_arr2 real(dp), dimension(1:Nz):: gxx_out,gxy_out,gxz_out,gyy_out,gyz_out,gzz_out,Bfield_out,jacobian_out, dBdx_out, dBdz_out, chi_out real(dp), dimension(1:Nz):: R_out, Z_out, dxdR_out, dxdZ_out real(dp):: d_inv, drPsi, dxPsi, dq_dx, dq_dpsi, R0, Z0, B0, F, D0_full, D1_full, D2_full, D3_full !real(dp) :: Lnorm, Psi0 ! currently module-wide defined anyway real(dp):: pprime, ffprime, D0_mid, D1_mid, D2_mid, D3_mid, dx_drho, pi, mu_0, dzprimedz real(dp):: rho_a, psiN, drpsiN, CN2, CN3, Rcenter, Zcenter, drRcenter, drZcenter logical:: bMaxShift integer:: i, k, iBmax Npol_ext = Npol+2 Npol_s = Npol+1 np = 500*Npol_ext np_s = 500*Npol_s rho = trpeps*major_R if (rho.le.0.0) stop 'flux surface radius not defined' trpeps = rho/major_R q0 = sign_Ip_CW * sign_Bt_CW * abs(q0) R0=major_R B0=1.0*sign_Bt_CW F=R0*B0 Z0=major_Z pi = acos(-1.0) mu_0=4.0*pi theta=linspace(-pi*Npol_ext,pi*Npol_ext-2._dp*pi*Npol_ext/np,np) d_inv=asin(delta) thetaShift = 0.0 iBmax = 1 !flux surface parametrization thAdj = theta + thetaShift if (zeta/=0.0 .or. s_zeta/=0.0) then R = R0 + rho*Cos(thAdj + d_inv*Sin(thAdj)) Z = Z0 + kappa*rho*Sin(thAdj + zeta*Sin(2*thAdj)) R_rho = drR + Cos(thAdj + d_inv*Sin(thAdj)) - s_delta*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj)) Z_rho = drZ + kappa*s_zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) & + kappa*Sin(thAdj + zeta*Sin(2*thAdj)) + kappa*s_kappa*Sin(thAdj + zeta*Sin(2*thAdj)) R_theta = -(rho*(1 + d_inv*Cos(thAdj))*Sin(thAdj + d_inv*Sin(thAdj))) Z_theta = kappa*rho*(1 + 2*zeta*Cos(2*thAdj))*Cos(thAdj + zeta*Sin(2*thAdj)) R_theta_theta = -(rho*(1 + d_inv*Cos(thAdj))**2*Cos(thAdj + d_inv*Sin(thAdj))) & + d_inv*rho*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj)) Z_theta_theta = -4*kappa*rho*zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) & - kappa*rho*(1 + 2*zeta*Cos(2*thAdj))**2*Sin(thAdj + zeta*Sin(2*thAdj)) else Rcirc = rho*Cos(thAdj - thetad + thetak) Zcirc = rho*Sin(thAdj - thetad + thetak) Relong = Rcirc Zelong = Zcirc + (-1 + kappa)*rho*Sin(thAdj - thetad + thetak) RelongTilt = Relong*Cos(thetad - thetak) - Zelong*Sin(thetad - thetak) ZelongTilt = Zelong*Cos(thetad - thetak) + Relong*Sin(thetad - thetak) Rtri = RelongTilt - rho*Cos(thAdj) + rho*Cos(thAdj + delta*Sin(thAdj)) Ztri = ZelongTilt RtriTilt = Rtri*Cos(thetad) + Ztri*Sin(thetad) ZtriTilt = Ztri*Cos(thetad) - Rtri*Sin(thetad) R = R0 + RtriTilt Z = Z0 + ZtriTilt drRcirc = Cos(thAdj - thetad + thetak) drZcirc = Sin(thAdj - thetad + thetak) drRelong = drRcirc drZelong = drZcirc - (1 - kappa - kappa*s_kappa)*Sin(thAdj - thetad + thetak) drRelongTilt = drRelong*Cos(thetad - thetak) - drZelong*Sin(thetad - thetak) drZelongTilt = drZelong*Cos(thetad - thetak) + drRelong*Sin(thetad - thetak) drRtri = drRelongTilt - Cos(thAdj) + Cos(thAdj + delta*Sin(thAdj)) & - s_delta*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj)) drZtri = drZelongTilt drRtriTilt = drRtri*Cos(thetad) + drZtri*Sin(thetad) drZtriTilt = drZtri*Cos(thetad) - drRtri*Sin(thetad) R_rho = drR + drRtriTilt Z_rho = drZ + drZtriTilt dtRcirc = -(rho*Sin(thAdj - thetad + thetak)) dtZcirc = rho*Cos(thAdj - thetad + thetak) dtRelong = dtRcirc dtZelong = dtZcirc + (-1 + kappa)*rho*Cos(thAdj - thetad + thetak) dtRelongTilt = dtRelong*Cos(thetad - thetak) - dtZelong*Sin(thetad - thetak) dtZelongTilt = dtZelong*Cos(thetad - thetak) + dtRelong*Sin(thetad - thetak) dtRtri = dtRelongTilt + rho*Sin(thAdj) - rho*(1 + delta*Cos(thAdj))*Sin(thAdj + delta*Sin(thAdj)) dtZtri = dtZelongTilt dtRtriTilt = dtRtri*Cos(thetad) + dtZtri*Sin(thetad) dtZtriTilt = dtZtri*Cos(thetad) - dtRtri*Sin(thetad) R_theta = dtRtriTilt Z_theta = dtZtriTilt dtdtRcirc = -(rho*Cos(thAdj - thetad + thetak)) dtdtZcirc = -(rho*Sin(thAdj - thetad + thetak)) dtdtRelong = dtdtRcirc dtdtZelong = dtdtZcirc - (-1 + kappa)*rho*Sin(thAdj - thetad + thetak) dtdtRelongTilt = dtdtRelong*Cos(thetad - thetak) - dtdtZelong*Sin(thetad - thetak) dtdtZelongTilt = dtdtZelong*Cos(thetad - thetak) + dtdtRelong*Sin(thetad - thetak) dtdtRtri = dtdtRelongTilt + rho*Cos(thAdj) - rho*(1 + delta*Cos(thAdj))**2*Cos(thAdj + delta*Sin(thAdj)) & + delta*rho*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj)) dtdtZtri = dtdtZelongTilt dtdtRtriTilt = dtdtRtri*Cos(thetad) + dtdtZtri*Sin(thetad) dtdtZtriTilt = dtdtZtri*Cos(thetad) - dtdtRtri*Sin(thetad) R_theta_theta = dtdtRtriTilt Z_theta_theta = dtdtZtriTilt endif !dl/dtheta dlp=(R_theta**2+Z_theta**2)**0.5 !curvature radius Rc=dlp**3*(R_theta*Z_theta_theta-Z_theta*R_theta_theta)**(-1) ! some useful quantities (see papers for definition of u) cosu=Z_theta/dlp sinu=-R_theta/dlp !Jacobian J_r = (dPsi/dr) J_psi = (dPsi/dr) / [(nabla fz x nabla psi)* nabla theta] ! = R * (dR/drho dZ/dtheta - dR/dtheta dZ/drho) = R dlp / |nabla r| J_r=R*(R_rho*Z_theta-R_theta*Z_rho) !From definition of q = 1/(2 pi) int (B nabla fz) / (B nabla theta) dtheta: !dPsi/dr = sign_Bt sign_Ip / (2 pi q) int F / R^2 J_r dtheta ! = F / (2 pi |q|) int J_r/R^2 dtheta tmp_arr=J_r/R**2 drPsi=sign_Ip_CW*F/(2.*pi*Npol_ext*q0)*sum(tmp_arr)*2*pi*Npol_ext/np !dlp_int(tmp_arr,1.0) !Poloidal field (Bp = Bvec * nabla l) Bp=sign_Ip_CW * drPsi / J_r * dlp !toroidal field Bphi=F/R !total modulus of Bfield B=sqrt(Bphi**2+Bp**2) bMaxShift = .false. ! if (thetaShift==0.0.and.trim(magn_geometry).ne.'miller_general') then if (thetaShift==0.0) then do i = 2,np-1 if (B(iBmax) x = r dx_drho=1. !drPsi/Psi0*Lnorm*q0 if (my_id==0) write(*,"(A,ES12.4)") 'Using radial coordinate with dx/dr = ',dx_drho dxPsi = drPsi/dx_drho C_y = dxPsi*sign_Ip_CW C_xy = abs(B0*dxPsi/C_y) if (my_id==0) then write(*,"(A,ES12.4,A,ES12.4,A,ES12.4)") & "Setting C_xy = ",C_xy,' C_y = ', C_y," C_x' = ", 1./dxPsi write(*,'(A,ES12.4)') "B_unit/Bref conversion factor = ", q0/rho*drPsi write(*,'(A,ES12.4)') "dPsi/dr = ", drPsi if (thetaShift.ne.0.0) write(*,'(A,ES12.4)') "thetaShift = ", thetaShift endif !--------shear is expected to be defined as rho/q*dq/drho--------! dq_dx=shat*q0/rho/dx_drho dq_dpsi=dq_dx/dxPsi pprime=-amhd/q0**2/R0/(2*mu_0)*B0**2/drPsi !neg. dpdx normalized to magnetic pressure for pressure term dpdx_pm_geom=amhd/q0**2/R0/dx_drho !first coefficient of psi in varrho expansion psi1 = R*Bp*sign_Ip_CW !integrals for ffprime evaluation do i=1,np tmp_arr=(2./Rc-2.*cosu/R)/(R*psi1**2) D0(i)=-F*dlp_int_ind(tmp_arr,dlp,i) tmp_arr=B**2*R/psi1**3 D1(i)=-dlp_int_ind(tmp_arr,dlp,i)/F tmp_arr=mu_0*R/psi1**3 D2(i)=-dlp_int_ind(tmp_arr,dlp,i)*F tmp_arr=1./(R*psi1) D3(i)=-dlp_int_ind(tmp_arr,dlp,i)*F enddo tmp_arr=(2./Rc-2.*cosu/R)/(R*psi1**2) D0_full=-F*dlp_int(tmp_arr,dlp) tmp_arr=B**2*R/psi1**3 D1_full=-dlp_int(tmp_arr,dlp)/F tmp_arr=mu_0*R/psi1**3 D2_full=-dlp_int(tmp_arr,dlp)*F tmp_arr=1./(R*psi1) D3_full=-dlp_int(tmp_arr,dlp)*F D0_mid=D0(np/2+1) D1_mid=D1(np/2+1) D2_mid=D2(np/2+1) D3_mid=D3(np/2+1) ffprime=-(sign_Ip_CW*dq_dpsi*2.*pi*Npol_ext+D0_full+D2_full*pprime)/D1_full if (my_id==0) then write(*,'(A,ES12.4)') "ffprime = ", ffprime endif D0=D0-D0_mid D1=D1-D1_mid D2=D2-D2_mid nu=D3-D3_mid nu1=psi1*(D0+D1*ffprime+D2*pprime) !straight field line angle defined on equidistant theta grid !alpha = fz + nu = - (q chi - fz) => chi = -nu / q chi=-nu/q0 !correct small scaling error (<0.5%, due to finite integration resolution) chi=chi*(maxval(theta)-minval(theta))/(maxval(chi)-minval(chi)) !new grid equidistant in straight field line angle chi_s = linspace(-pi*Npol_s,pi*Npol_s-2*pi*Npol_s/np_s,np_s) if (sign_Ip_CW.lt.0.0) then !make chi increasing function to not confuse lag3interp tmp_reverse = chi(np:1:-1) theta_reverse = theta(np:1:-1) call lag3interp(theta_reverse,tmp_reverse,np,theta_s,chi_s,np_s) theta_s_reverse = theta_s(np_s:1:-1) else !lag3interp(y_in,x_in,n_in,y_out,x_out,n_out) call lag3interp(theta,chi,np,theta_s,chi_s,np_s) endif dtheta_dchi_s=deriv_fd(theta_s,chi_s,np_s) !arrays equidistant in straight field line angle thAdj_s = theta_s + thetaShift if (zeta/=0.0 .or. s_zeta/=0.0) then R_s = R0 + rho*Cos(thAdj_s + d_inv*Sin(thAdj_s)) Z_s = Z0 + kappa*rho*Sin(thAdj_s + zeta*Sin(2*thAdj_s)) R_theta_s = -(dtheta_dchi_s*rho*(1 + d_inv*Cos(thAdj_s))*Sin(thAdj_s + d_inv*Sin(thAdj_s))) Z_theta_s = dtheta_dchi_s*kappa*rho*(1 + 2*zeta*Cos(2*thAdj_s))*Cos(thAdj_s + zeta*Sin(2*thAdj_s)) else Rcirc_s = rho*Cos(thAdj_s - thetad + thetak) Zcirc_s = rho*Sin(thAdj_s - thetad + thetak) Relong_s = Rcirc_s Zelong_s = Zcirc_s + (-1 + kappa)*rho*Sin(thAdj_s - thetad + thetak) RelongTilt_s = Relong_s*Cos(thetad - thetak) - Zelong_s*Sin(thetad - thetak) ZelongTilt_s = Zelong_s*Cos(thetad - thetak) + Relong_s*Sin(thetad - thetak) Rtri_s = RelongTilt_s - rho*Cos(thAdj_s) + rho*Cos(thAdj_s + delta*Sin(thAdj_s)) Ztri_s = ZelongTilt_s RtriTilt_s = Rtri_s*Cos(thetad) + Ztri_s*Sin(thetad) ZtriTilt_s = Ztri_s*Cos(thetad) - Rtri_s*Sin(thetad) R_s = R0 + RtriTilt_s Z_s = Z0 + ZtriTilt_s dtRcirc_s = -(rho*Sin(thAdj_s - thetad + thetak)) dtZcirc_s = rho*Cos(thAdj_s - thetad + thetak) dtRelong_s = dtRcirc_s dtZelong_s = dtZcirc_s + (-1 + kappa)*rho*Cos(thAdj_s - thetad + thetak) dtRelongTilt_s = dtRelong_s*Cos(thetad - thetak) - dtZelong_s*Sin(thetad - thetak) dtZelongTilt_s = dtZelong_s*Cos(thetad - thetak) + dtRelong_s*Sin(thetad - thetak) dtRtri_s = dtRelongTilt_s + rho*Sin(thAdj_s) & - rho*(1 + delta*Cos(thAdj_s))*Sin(thAdj_s + delta*Sin(thAdj_s)) dtZtri_s = dtZelongTilt_s dtRtriTilt_s = dtRtri_s*Cos(thetad) + dtZtri_s*Sin(thetad) dtZtriTilt_s = dtZtri_s*Cos(thetad) - dtRtri_s*Sin(thetad) R_theta_s = dtheta_dchi_s*dtRtriTilt_s Z_theta_s = dtheta_dchi_s*dtZtriTilt_s endif if (sign_Ip_CW.lt.0.0) then call lag3interp(nu1,theta,np,tmp_arr_s,theta_s_reverse,np_s) nu1_s = tmp_arr_s(np_s:1:-1) call lag3interp(Bp,theta,np,tmp_arr_s,theta_s_reverse,np_s) Bp_s = tmp_arr_s(np_s:1:-1) call lag3interp(dlp,theta,np,tmp_arr_s,theta_s_reverse,np_s) dlp_s = tmp_arr_s(np_s:1:-1) call lag3interp(Rc,theta,np,tmp_arr_s,theta_s_reverse,np_s) Rc_s = tmp_arr_s(np_s:1:-1) else call lag3interp(nu1,theta,np,nu1_s,theta_s,np_s) call lag3interp(Bp,theta,np,Bp_s,theta_s,np_s) call lag3interp(dlp,theta,np,dlp_s,theta_s,np_s) call lag3interp(Rc,theta,np,Rc_s,theta_s,np_s) endif psi1_s = R_s*Bp_s*sign_Ip_CW dBp_dchi_s=deriv_fd(Bp_s,chi_s,np_s) Bphi_s=F/R_s B_s=sqrt(Bphi_s**2+Bp_s**2) cosu_s=Z_theta_s/dlp_s/dtheta_dchi_s sinu_s=-R_theta_s/dlp_s/dtheta_dchi_s !radial derivative of straight field line angle dchidx_s=-(nu1_s/psi1_s*dxPsi+chi_s*dq_dx)/q0 !Bfield derivatives in Mercier-Luc coordinates (varrho,l,fz) dB_drho_s=-1./B_s*(F**2/R_s**3*cosu_s+Bp_s**2/Rc_s+mu_0*psi1_s*pprime) dB_dl_s=1./B_s*(Bp_s*dBp_dchi_s/dtheta_dchi_s/dlp_s+F**2/R_s**3*sinu_s) dnu_drho_s=nu1_s dnu_dl_s=-F/(R_s*psi1_s) grad_nu_s=sqrt(dnu_drho_s**2+dnu_dl_s**2) !contravariant metric coefficients (varrho,l,fz)->(x,y,z) gxx=(psi1_s/dxPsi)**2 gxy=-psi1_s/dxPsi*C_y*sign_Ip_CW*nu1_s gxz=-psi1_s/dxPsi*(nu1_s+psi1_s*dq_dpsi*chi_s)/q0 gyy=C_y**2*(grad_nu_s**2+1/R_s**2) gyz=sign_Ip_CW*C_y/q0*(grad_nu_s**2+dq_dpsi*nu1_s*psi1_s*chi_s) gzz=1./q0**2*(grad_nu_s**2+2.*dq_dpsi*nu1_s*psi1_s*chi_s+(dq_dpsi*psi1_s*chi_s)**2) jacobian=1./sqrt(gxx*gyy*gzz + 2.*gxy*gyz*gxz - gxz**2*gyy - gyz**2*gxx - gzz*gxy**2) !covariant metric coefficients g_xx=jacobian**2*(gyy*gzz-gyz**2) g_xy=jacobian**2*(gxz*gyz-gxy*gzz) g_xz=jacobian**2*(gxy*gyz-gxz*gyy) g_yy=jacobian**2*(gxx*gzz-gxz**2) g_yz=jacobian**2*(gxz*gxy-gxx*gyz) g_zz=jacobian**2*(gxx*gyy-gxy**2) !Bfield derivatives !dBdx = e_x * nabla B = J (nabla y x nabla z) * nabla B dBdx=jacobian*C_y/(q0*R_s)*(F/(R_s*psi1_s)*dB_drho_s+(nu1_s+dq_dpsi*chi_s*psi1_s)*dB_dl_s) dBdz=1./B_s*(Bp_s*dBp_dchi_s-F**2/R_s**3*R_theta_s) !curvature terms (these are just local and will be recalculated in geometry.F90) K_x = (0.-g_yz/g_zz*dBdz) K_y = (dBdx-g_xz/g_zz*dBdz) !(R,Z) derivatives for visualization dxdR_s = dx_drho/drPsi*psi1_s*cosu_s dxdZ_s = dx_drho/drPsi*psi1_s*sinu_s if (edge_opt==0.0) then !gene z-grid chi_out=linspace(-pi*Npol,pi*Npol-2*pi*Npol/Nz,Nz) else !new parallel coordinate chi_out==zprime !see also tracer_aux.F90 if (Npol>1) STOP "ERROR: Npol>1 has not been implemented for edge_opt=\=0.0" do k=izs,ize chi_out(k)=sinh((-pi+k*2.*pi/Nz)*log(edge_opt*pi+sqrt(edge_opt**2*pi**2+1))/pi)/edge_opt enddo !transform metrics according to chain rule do k=1,np_s !>dz'/dz conversion for edge_opt as function of z if (edge_opt.gt.0) then dzprimedz = edge_opt*pi/log(edge_opt*pi+sqrt((edge_opt*pi)**2+1))/& sqrt((edge_opt*chi_s(k))**2+1) else dzprimedz = 1.0 endif gxz(k)=gxz(k)*dzprimedz gyz(k)=gyz(k)*dzprimedz gzz(k)=gzz(k)*dzprimedz**2 jacobian(k)=jacobian(k)/dzprimedz dBdz(k)=dBdz(k)/dzprimedz enddo endif !edge_opt !interpolate down to GENE z-grid call lag3interp(gxx,chi_s,np_s,gxx_out,chi_out,Nz) call lag3interp(gxy,chi_s,np_s,gxy_out,chi_out,Nz) call lag3interp(gxz,chi_s,np_s,gxz_out,chi_out,Nz) call lag3interp(gyy,chi_s,np_s,gyy_out,chi_out,Nz) call lag3interp(gyz,chi_s,np_s,gyz_out,chi_out,Nz) call lag3interp(gzz,chi_s,np_s,gzz_out,chi_out,Nz) call lag3interp(B_s,chi_s,np_s,Bfield_out,chi_out,Nz) call lag3interp(jacobian,chi_s,np_s,jacobian_out,chi_out,Nz) call lag3interp(dBdx,chi_s,np_s,dBdx_out,chi_out,Nz) call lag3interp(dBdz,chi_s,np_s,dBdz_out,chi_out,Nz) call lag3interp(R_s,chi_s,np_s,R_out,chi_out,Nz) call lag3interp(Z_s,chi_s,np_s,Z_out,chi_out,Nz) call lag3interp(dxdR_s,chi_s,np_s,dxdR_out,chi_out,Nz) call lag3interp(dxdZ_s,chi_s,np_s,dxdZ_out,chi_out,Nz) ! Fill the geom arrays with the results do eo=0,1 gxx_(izs:ize,eo) =gxx_out(izs:ize) gyy_(izs:ize,eo) =gyy_out(izs:ize) gxz_(izs:ize,eo) =gxz_out(izs:ize) gyz_(izs:ize,eo) =gyz_out(izs:ize) dBdx_(izs:ize,eo) =dBdx_out(izs:ize) dBdy_(izs:ize,eo) =0. gxy_(izs:ize,eo) =gxy_out(izs:ize) gzz_(izs:ize,eo) =gzz_out(izs:ize) Bfield_(izs:ize,eo) =Bfield_out(izs:ize) jacobian_(izs:ize,eo) =jacobian_out(izs:ize) dBdz_(izs:ize,eo) =dBdz_out(izs:ize) R_hat_(izs:ize,eo) =R_out(izs:ize) Z_hat_(izs:ize,eo) =Z_out(izs:ize) dxdR_(izs:ize,eo) = dxdR_out(izs:ize) dxdZ_(izs:ize,eo) = dxdZ_out(izs:ize) !! UPDATE GHOSTS VALUES (since the miller function in GENE does not) CALL update_ghosts_z(gxx_(:,eo)) CALL update_ghosts_z(gyy_(:,eo)) CALL update_ghosts_z(gxz_(:,eo)) - CALL update_ghosts_z(dBdx_(:,eo)) - CALL update_ghosts_z(dBdy_(:,eo)) CALL update_ghosts_z(gxy_(:,eo)) CALL update_ghosts_z(gzz_(:,eo)) CALL update_ghosts_z(Bfield_(:,eo)) - CALL update_ghosts_z(jacobian_(:,eo)) + CALL update_ghosts_z(dBdx_(:,eo)) + CALL update_ghosts_z(dBdy_(:,eo)) CALL update_ghosts_z(dBdz_(:,eo)) + CALL update_ghosts_z(jacobian_(:,eo)) CALL update_ghosts_z(R_hat_(:,eo)) CALL update_ghosts_z(Z_hat_(:,eo)) CALL update_ghosts_z(dxdR_(:,eo)) CALL update_ghosts_z(dxdZ_(:,eo)) + + ! Curvature operator (Frei et al. 2022 eq 2.15) + DO iz = izgs,izge + G1 = gxx_(iz,eo)*gyy_(iz,eo)-gxy_(iz,eo)*gxy_(iz,eo) + G2 = gxx_(iz,eo)*gyz_(iz,eo)-gxy_(iz,eo)*gxz_(iz,eo) + G3 = gxy_(iz,eo)*gyz_(iz,eo)-gyy_(iz,eo)*gxz_(iz,eo) + Cx = (G1*dBdy_(iz,eo) + G2*dBdz_(iz,eo))/Bfield_(iz,eo) + Cy = (G3*dBdz_(iz,eo) - G1*dBdx_(iz,eo))/Bfield_(iz,eo) + + DO iky = ikys, ikye + ky = kyarray(iky) + DO ikx= ikxs, ikxe + kx = kxarray(ikx) + Ckxky_(iky, ikx, iz,eo) = (Cx*kx + Cy*ky)/Bfield_(iz,eo) + ENDDO + ENDDO + ENDDO + ! coefficient in the front of parallel derivative + gradz_coeff_(iz,eo) = 1._dp / Jacobian_(iz,eo) / Bfield_(iz,eo) enddo contains SUBROUTINE update_ghosts_z(fz_) IMPLICIT NONE ! INTEGER, INTENT(IN) :: nztot_ REAL(dp), DIMENSION(izgs:izge), INTENT(INOUT) :: fz_ REAL(dp), DIMENSION(-2:2) :: buff INTEGER :: status(MPI_STATUS_SIZE), count IF(Nz .GT. 1) THEN IF (num_procs_z .GT. 1) THEN CALL MPI_BARRIER(MPI_COMM_WORLD,ierr) count = 1 ! one point to exchange !!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!! CALL mpi_sendrecv(fz_(ize), count, MPI_DOUBLE, nbr_U, 0, & ! Send to Up the last buff(-1), count, MPI_DOUBLE, nbr_D, 0, & ! Receive from Down the first-1 comm0, status, ierr) CALL mpi_sendrecv(fz_(ize-1), count, MPI_DOUBLE, nbr_U, 0, & ! Send to Up the last buff(-2), count, MPI_DOUBLE, nbr_D, 0, & ! Receive from Down the first-2 comm0, status, ierr) !!!!!!!!!!! Send ghost to down neighbour !!!!!!!!!!!!!!!!!!!!!! CALL mpi_sendrecv(fz_(izs), count, MPI_DOUBLE, nbr_D, 0, & ! Send to Down the first buff(+1), count, MPI_DOUBLE, nbr_U, 0, & ! Recieve from Up the last+1 comm0, status, ierr) CALL mpi_sendrecv(fz_(izs+1), count, MPI_DOUBLE, nbr_D, 0, & ! Send to Down the first buff(+2), count, MPI_DOUBLE, nbr_U, 0, & ! Recieve from Up the last+2 comm0, status, ierr) ELSE buff(-1) = fz_(ize ) buff(-2) = fz_(ize-1) buff(+1) = fz_(izs ) buff(+2) = fz_(izs+1) ENDIF fz_(ize+1) = buff(+1) fz_(ize+2) = buff(+2) fz_(izs-1) = buff(-1) fz_(izs-2) = buff(-2) ENDIF END SUBROUTINE update_ghosts_z !> Generate an equidistant array from min to max with n points function linspace(min,max,n) result(out) real(dp):: min, max integer:: n real(dp), dimension(n):: out do i=1,n out(i)=min+(i-1)*(max-min)/(n-1) enddo end function linspace !> Weighted average real(dp) function average(var,weight) real(dp), dimension(np):: var, weight average=sum(var*weight)/sum(weight) end function average !> full theta integral with weight function dlp real(dp) function dlp_int(var,dlp) real(dp), dimension(np):: var, dlp dlp_int=sum(var*dlp)*2*pi*Npol_ext/np end function dlp_int !> theta integral with weight function dlp, up to index 'ind' real(dp) function dlp_int_ind(var,dlp,ind) real(dp), dimension(np):: var, dlp integer:: ind dlp_int_ind=0. if (ind.gt.1) then dlp_int_ind=dlp_int_ind+var(1)*dlp(1)*pi*Npol_ext/np dlp_int_ind=dlp_int_ind+(sum(var(2:ind-1)*dlp(2:ind-1)))*2*pi*Npol_ext/np dlp_int_ind=dlp_int_ind+var(ind)*dlp(ind)*pi*Npol_ext/np endif end function dlp_int_ind !> 1st derivative with 2nd order finite differences function deriv_fd(y,x,n) result(out) integer, intent(in) :: n real(dp), dimension(n):: x,y,out,dx !call lag3deriv(y,x,n,out,x,n) out=0. do i=2,n-1 out(i)=out(i)-y(i-1)/2 out(i)=out(i)+y(i+1)/2 enddo out(1)=y(2)-y(1) out(n)=y(n)-y(n-1) dx=x(2)-x(1) out=out/dx end function deriv_fd end subroutine get_miller END MODULE miller diff --git a/src/model_mod.F90 b/src/model_mod.F90 index 75a2768..5d28283 100644 --- a/src/model_mod.F90 +++ b/src/model_mod.F90 @@ -1,152 +1,154 @@ MODULE model ! Module for diagnostic parameters USE prec_const IMPLICIT NONE PRIVATE INTEGER, PUBLIC, PROTECTED :: CLOS = 0 ! linear truncation method INTEGER, PUBLIC, PROTECTED :: NL_CLOS = 0 ! nonlinear truncation method INTEGER, PUBLIC, PROTECTED :: KERN = 0 ! Kernel model CHARACTER(len=16), & PUBLIC, PROTECTED ::LINEARITY= 'linear' ! To turn on non linear bracket term LOGICAL, PUBLIC, PROTECTED :: KIN_E = .true. ! Simulate kinetic electron (adiabatic otherwise) REAL(dp), PUBLIC, PROTECTED :: mu_x = 0._dp ! spatial x-Hyperdiffusivity coefficient (for num. stability) REAL(dp), PUBLIC, PROTECTED :: mu_y = 0._dp ! spatial y-Hyperdiffusivity coefficient (for num. stability) REAL(dp), PUBLIC, PROTECTED :: mu_z = 0._dp ! spatial z-Hyperdiffusivity coefficient (for num. stability) REAL(dp), PUBLIC, PROTECTED :: mu_p = 0._dp ! kinetic para hyperdiffusivity coefficient (for num. stability) REAL(dp), PUBLIC, PROTECTED :: mu_j = 0._dp ! kinetic perp hyperdiffusivity coefficient (for num. stability) INTEGER, PUBLIC, PROTECTED :: N_HD = 4 ! order of numerical spatial diffusion REAL(dp), PUBLIC, PROTECTED :: nu = 0._dp ! Collision frequency REAL(dp), PUBLIC, PROTECTED :: tau_e = 1._dp ! Temperature REAL(dp), PUBLIC, PROTECTED :: tau_i = 1._dp ! REAL(dp), PUBLIC, PROTECTED :: sigma_e = 0.023338_dp! sqrt of electron ion mass ratio REAL(dp), PUBLIC, PROTECTED :: sigma_i = 1._dp ! REAL(dp), PUBLIC, PROTECTED :: q_e = -1._dp ! Charge REAL(dp), PUBLIC, PROTECTED :: q_i = 1._dp ! REAL(dp), PUBLIC, PROTECTED :: k_Ni = 7._dp ! Ion density drive REAL(dp), PUBLIC, PROTECTED :: k_Ne = 7._dp ! Ele '' REAL(dp), PUBLIC, PROTECTED :: k_Ti = 2._dp ! Ion temperature drive REAL(dp), PUBLIC, PROTECTED :: k_Te = 2._dp ! Ele '' REAL(dp), PUBLIC, PROTECTED :: K_E = 0._dp ! Backg. electric field drive REAL(dp), PUBLIC, PROTECTED :: GradB = 1._dp ! Magnetic gradient REAL(dp), PUBLIC, PROTECTED :: CurvB = 1._dp ! Magnetic curvature REAL(dp), PUBLIC, PROTECTED :: lambdaD = 1._dp ! Debye length REAL(dp), PUBLIC, PROTECTED :: beta = 0._dp ! electron plasma Beta (8piNT_e/B0^2) LOGICAL, PUBLIC, PROTECTED :: EM = .false. ! Electromagnetic effects flag REAL(dp), PUBLIC, PROTECTED :: taue_qe ! factor of the magnetic moment coupling REAL(dp), PUBLIC, PROTECTED :: taui_qi ! REAL(dp), PUBLIC, PROTECTED :: qi_taui ! REAL(dp), PUBLIC, PROTECTED :: qe_taue ! REAL(dp), PUBLIC, PROTECTED :: sqrtTaue_qe ! factor of parallel moment term REAL(dp), PUBLIC, PROTECTED :: sqrtTaui_qi ! REAL(dp), PUBLIC, PROTECTED :: qe_sigmae_sqrtTaue ! factor of parallel phi term REAL(dp), PUBLIC, PROTECTED :: qi_sigmai_sqrtTaui ! REAL(dp), PUBLIC, PROTECTED :: sigmae2_taue_o2 ! factor of the Kernel argument REAL(dp), PUBLIC, PROTECTED :: sigmai2_taui_o2 ! REAL(dp), PUBLIC, PROTECTED :: sqrt_sigmae2_taue_o2 ! factor of the Kernel argument REAL(dp), PUBLIC, PROTECTED :: sqrt_sigmai2_taui_o2 REAL(dp), PUBLIC, PROTECTED :: nu_e, nu_i ! electron-ion, ion-ion collision frequency REAL(dp), PUBLIC, PROTECTED :: nu_ee, nu_ie ! e-e, i-e coll. frequ. REAL(dp), PUBLIC, PROTECTED :: qe2_taue, qi2_taui ! factor of the gammaD sum REAL(dp), PUBLIC, PROTECTED :: q_o_sqrt_tau_sigma_e, q_o_sqrt_tau_sigma_i REAL(dp), PUBLIC, PROTECTED :: sqrt_tau_o_sigma_e, sqrt_tau_o_sigma_i + REAL(dp), PUBLIC, PROTECTED :: dpdx = 0 ! radial pressure gradient PUBLIC :: model_readinputs, model_outputinputs CONTAINS SUBROUTINE model_readinputs ! Read the input parameters USE basic, ONLY : lu_in, my_id USE prec_const IMPLICIT NONE NAMELIST /MODEL_PAR/ CLOS, NL_CLOS, KERN, LINEARITY, KIN_E, & mu_x, mu_y, N_HD, mu_z, mu_p, mu_j, nu,& tau_e, tau_i, sigma_e, sigma_i, q_e, q_i,& k_Ne, k_Ni, k_Te, k_Ti, GradB, CurvB, lambdaD, beta READ(lu_in,model_par) IF(.NOT. KIN_E) THEN IF(my_id.EQ.0) print*, 'Adiabatic electron model -> beta = 0' beta = 0._dp ENDIF IF(beta .GT. 0) THEN - EM = .TRUE. IF(my_id.EQ.0) print*, 'Electromagnetic effects are included' + EM = .TRUE. + dpdx = -(tau_i*(k_Ni + k_Ti) + tau_e*(k_Ne + k_Te)) ENDIF taue_qe = tau_e/q_e ! factor of the magnetic moment coupling taui_qi = tau_i/q_i ! factor of the magnetic moment coupling qe_taue = q_e/tau_e qi_taui = q_i/tau_i sqrtTaue_qe = sqrt(tau_e)/q_e ! factor of parallel moment term sqrtTaui_qi = sqrt(tau_i)/q_i ! factor of parallel moment term qe_sigmae_sqrtTaue = q_e/sigma_e/SQRT(tau_e) ! factor of parallel phi term qi_sigmai_sqrtTaui = q_i/sigma_i/SQRT(tau_i) qe2_taue = (q_e**2)/tau_e ! factor of the gammaD sum qi2_taui = (q_i**2)/tau_i sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp sqrt_sigmae2_taue_o2 = SQRT(sigma_e**2 * tau_e/2._dp) ! to avoid multiple SQRT eval sqrt_sigmai2_taui_o2 = SQRT(sigma_i**2 * tau_i/2._dp) q_o_sqrt_tau_sigma_e = q_e/SQRT(tau_e)/sigma_e ! For psi field terms q_o_sqrt_tau_sigma_i = q_i/SQRT(tau_i)/sigma_i ! For psi field terms sqrt_tau_o_sigma_e = SQRT(tau_e)/sigma_e ! For Ampere eq sqrt_tau_o_sigma_i = SQRT(tau_i)/sigma_i !! We use the ion-ion collision as normalization with definition ! nu_ii = 4 sqrt(pi)/3 T_i^(-3/2) m_i^(-1/2) q^4 n_i0 ln(Lambda) ! nu_e = nu/sigma_e * (tau_e)**(3._dp/2._dp) ! electron-ion collision frequency (where already multiplied by 0.532) nu_i = nu ! ion-ion collision frequ. nu_ee = nu_e ! e-e coll. frequ. nu_ie = nu_i ! i-e coll. frequ. ! Old normalization (MOLI Jorge/Frei) ! nu_e = 0.532_dp*nu ! electron-ion collision frequency (where already multiplied by 0.532) ! nu_i = 0.532_dp*nu*sigma_e*tau_e**(-3._dp/2._dp)/SQRT2 ! ion-ion collision frequ. ! nu_ee = nu_e/SQRT2 ! e-e coll. frequ. ! nu_ie = 0.532_dp*nu*sigma_e**2 ! i-e coll. frequ. END SUBROUTINE model_readinputs SUBROUTINE model_outputinputs(fidres, str) ! Write the input parameters to the results_xx.h5 file USE futils, ONLY: attach USE prec_const IMPLICIT NONE INTEGER, INTENT(in) :: fidres CHARACTER(len=256), INTENT(in) :: str CALL attach(fidres, TRIM(str), "CLOS", CLOS) CALL attach(fidres, TRIM(str), "KERN", KERN) CALL attach(fidres, TRIM(str), "LINEARITY", LINEARITY) CALL attach(fidres, TRIM(str), "KIN_E", KIN_E) CALL attach(fidres, TRIM(str), "nu", nu) CALL attach(fidres, TRIM(str), "mu", 0) CALL attach(fidres, TRIM(str), "mu_x", mu_x) CALL attach(fidres, TRIM(str), "mu_y", mu_y) CALL attach(fidres, TRIM(str), "mu_z", mu_z) CALL attach(fidres, TRIM(str), "mu_p", mu_p) CALL attach(fidres, TRIM(str), "mu_j", mu_j) CALL attach(fidres, TRIM(str), "tau_e", tau_e) CALL attach(fidres, TRIM(str), "tau_i", tau_i) CALL attach(fidres, TRIM(str), "sigma_e", sigma_e) CALL attach(fidres, TRIM(str), "sigma_i", sigma_i) CALL attach(fidres, TRIM(str), "q_e", q_e) CALL attach(fidres, TRIM(str), "q_i", q_i) CALL attach(fidres, TRIM(str), "k_Ne", k_Ne) CALL attach(fidres, TRIM(str), "k_Ni", k_Ni) CALL attach(fidres, TRIM(str), "k_Te", k_Te) CALL attach(fidres, TRIM(str), "k_Ti", k_Ti) CALL attach(fidres, TRIM(str), "K_E", K_E) CALL attach(fidres, TRIM(str), "lambdaD", lambdaD) CALL attach(fidres, TRIM(str), "beta", beta) END SUBROUTINE model_outputinputs END MODULE model diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90 index 8123fe7..948d66c 100644 --- a/src/moments_eq_rhs_mod.F90 +++ b/src/moments_eq_rhs_mod.F90 @@ -1,403 +1,407 @@ MODULE moments_eq_rhs IMPLICIT NONE PUBLIC :: compute_moments_eq_rhs CONTAINS SUBROUTINE compute_moments_eq_rhs USE model, only: KIN_E IMPLICIT NONE IF(KIN_E) CALL moments_eq_rhs_e CALL moments_eq_rhs_i !! BETA TESTING !! IF(.false.) CALL add_Maxwellian_background_terms END SUBROUTINE compute_moments_eq_rhs !_____________________________________________________________________________! !_____________________________________________________________________________! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! Electrons moments RHS !!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !_____________________________________________________________________________! SUBROUTINE moments_eq_rhs_e USE basic USE time_integration USE array, ONLY: xnepj, xnepp2j, xnepm2j, xnepjp1, xnepjm1, xnepp1j, xnepm1j,& ynepp1j, ynepp1jm1, ynepm1j, ynepm1jm1, & znepm1j, znepm1jp1, znepm1jm1, & xphij_e, xphijp1_e, xphijm1_e, xpsij_e, xpsijp1_e, xpsijm1_e,& nadiab_moments_e, ddz_nepj, interp_nepj, moments_rhs_e, Sepj,& TColl_e, ddzND_nepj, kernel_e USE fields, ONLY: phi, psi, moments_e USE grid USE model USE prec_const USE collision - USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatR + USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatB USE calculus, ONLY : interp_z, grad_z, grad_z2 IMPLICIT NONE INTEGER :: p_int, j_int ! loops indices and polynom. degrees REAL(dp) :: kx, ky, kperp2, dzlnB_o_J COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1 ! Terms from b x gradB and drives COMPLEX(dp) :: Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi, Tpsi COMPLEX(dp) :: Unepm1j, Unepm1jp1, Unepm1jm1 ! Terms from mirror force with adiab moments COMPLEX(dp) :: i_kx,i_ky,phikykxz, psikykxz ! Measuring execution time CALL cpu_time(t0_rhs) ! Spatial loops zloope : DO iz = izs,ize kxloope : DO ikx = ikxs,ikxe kx = kxarray(ikx) ! radial wavevector i_kx = imagu * kx ! toroidal derivative kyloope : DO iky = ikys,ikye ky = kyarray(iky) ! toroidal wavevector i_ky = imagu * ky ! toroidal derivative phikykxz = phi(iky,ikx,iz)! tmp phi value psikykxz = psi(iky,ikx,iz)! tmp psi value ! Kinetic loops jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxi+1 j_int = jarray_e(ij) ploope : DO ip = ips_e, ipe_e ! Hermite loop p_int = parray_e(ip) ! Hermite degree eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid kperp2= kparray(iky,ikx,iz,eo)**2 IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN !! Compute moments mixing terms Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp ! Perpendicular dynamic ! term propto n_e^{p,j} Tnepj = xnepj(ip,ij)* nadiab_moments_e(ip,ij,iky,ikx,iz) ! term propto n_e^{p+2,j} Tnepp2j = xnepp2j(ip) * nadiab_moments_e(ip+pp2,ij,iky,ikx,iz) ! term propto n_e^{p-2,j} Tnepm2j = xnepm2j(ip) * nadiab_moments_e(ip-pp2,ij,iky,ikx,iz) ! term propto n_e^{p,j+1} Tnepjp1 = xnepjp1(ij) * nadiab_moments_e(ip,ij+1,iky,ikx,iz) ! term propto n_e^{p,j-1} Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,iky,ikx,iz) ! Tperp Tperp = Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1 ! Parallel dynamic ! ddz derivative for Landau damping term Tpar = xnepp1j(ip) * ddz_nepj(ip+1,ij,iky,ikx,iz) & + xnepm1j(ip) * ddz_nepj(ip-1,ij,iky,ikx,iz) ! Mirror terms Tnepp1j = ynepp1j (ip,ij) * interp_nepj(ip+1,ij ,iky,ikx,iz) Tnepp1jm1 = ynepp1jm1(ip,ij) * interp_nepj(ip+1,ij-1,iky,ikx,iz) Tnepm1j = ynepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz) Tnepm1jm1 = ynepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz) ! Trapping terms Unepm1j = znepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz) Unepm1jp1 = znepm1jp1(ip,ij) * interp_nepj(ip-1,ij+1,iky,ikx,iz) Unepm1jm1 = znepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz) Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + Unepm1j + Unepm1jp1 + Unepm1jm1 !! Electrical potential term IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 Tphi = (xphij_e (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) & + xphijp1_e(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) & + xphijm1_e(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*phikykxz ELSE Tphi = 0._dp ENDIF !! Vector potential term IF ( (p_int .LE. 3) .AND. (p_int .GE. 1) ) THEN ! Kronecker p1 or p3 Tpsi = (xpsij_e (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) & + xpsijp1_e(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) & + xpsijm1_e(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*psikykxz ELSE Tpsi = 0._dp ENDIF !! Sum of all RHS terms moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = & ! Perpendicular magnetic gradient/curvature effects - -imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp& + -imagu*Ckxky(iky,ikx,iz,eo)*hatB(iz,eo) * Tperp& + ! Perpendicular pressure effects + -i_ky*beta*dpdx * Tperp& ! Parallel coupling (Landau Damping) -gradz_coeff(iz,eo) * Tpar & ! Mirror term (parallel magnetic gradient) -gradzB(iz,eo)*gradz_coeff(iz,eo) * Tmir& ! Drives (density + temperature gradients) -i_ky * (Tphi - Tpsi) & ! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose) -mu_x*diff_kx_coeff*kx**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) & -mu_y*diff_ky_coeff*ky**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) & ! Numerical parallel hyperdiffusion "mu_z*ddz**4" see Pueschel 2010 (eq 25) -mu_z*diff_dz_coeff*ddzND_nepj(ip,ij,iky,ikx,iz) & ! Collision term +TColl_e(ip,ij,iky,ikx,iz) & ! Nonlinear term -hatB_NL(iz,eo) * Sepj(ip,ij,iky,ikx,iz) IF(ip-4 .GT. 0) & ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33) moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = & moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) & + mu_p * moments_e(ip-4,ij,iky,ikx,iz,updatetlevel) ELSE moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp ENDIF END DO ploope END DO jloope END DO kyloope END DO kxloope END DO zloope ! Execution time end CALL cpu_time(t1_rhs) tc_rhs = tc_rhs + (t1_rhs-t0_rhs) END SUBROUTINE moments_eq_rhs_e !_____________________________________________________________________________! !_____________________________________________________________________________! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!! Ions moments RHS !!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !_____________________________________________________________________________! SUBROUTINE moments_eq_rhs_i USE basic USE time_integration, ONLY: updatetlevel USE array, ONLY: xnipj, xnipp2j, xnipm2j, xnipjp1, xnipjm1, xnipp1j, xnipm1j,& ynipp1j, ynipp1jm1, ynipm1j, ynipm1jm1, & znipm1j, znipm1jp1, znipm1jm1, & xphij_i, xphijp1_i, xphijm1_i, xpsij_i, xpsijp1_i, xpsijm1_i,& nadiab_moments_i, ddz_nipj, interp_nipj, moments_rhs_i, Sipj,& TColl_i, ddzND_nipj, kernel_i USE fields, ONLY: phi, psi, moments_i USE grid USE model USE prec_const USE collision - USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatR + USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatB USE calculus, ONLY : interp_z, grad_z, grad_z2 IMPLICIT NONE INTEGER :: p_int, j_int ! loops indices and polynom. degrees REAL(dp) :: kx, ky, kperp2 COMPLEX(dp) :: Tnipj, Tnipp2j, Tnipm2j, Tnipjp1, Tnipjm1 COMPLEX(dp) :: Tnipp1j, Tnipm1j, Tnipp1jm1, Tnipm1jm1 ! Terms from mirror force with non adiab moments COMPLEX(dp) :: Unipm1j, Unipm1jp1, Unipm1jm1 ! Terms from mirror force with adiab moments COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi, Tpsi COMPLEX(dp) :: i_kx, i_ky, phikykxz, psikykxz ! Measuring execution time CALL cpu_time(t0_rhs) ! Spatial loops zloopi : DO iz = izs,ize kxloopi : DO ikx = ikxs,ikxe kx = kxarray(ikx) ! radial wavevector i_kx = imagu * kx ! toroidal derivative kyloopi : DO iky = ikys,ikye ky = kyarray(iky) ! toroidal wavevector i_ky = imagu * ky ! toroidal derivative phikykxz = phi(iky,ikx,iz)! tmp phi value psikykxz = psi(iky,ikx,iz)! tmp phi value ! Kinetic loops jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1 j_int = jarray_i(ij) ploopi : DO ip = ips_i, ipe_i ! Hermite loop p_int = parray_i(ip) ! Hermite degree eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid kperp2= kparray(iky,ikx,iz,eo)**2 IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN !! Compute moments mixing terms Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp ! Perpendicular dynamic ! term propto n_i^{p,j} Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip ,ij ,iky,ikx,iz) ! term propto n_i^{p+2,j} Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij ,iky,ikx,iz) ! term propto n_i^{p-2,j} Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij ,iky,ikx,iz) ! term propto n_i^{p,j+1} Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip ,ij+1,iky,ikx,iz) ! term propto n_i^{p,j-1} Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip ,ij-1,iky,ikx,iz) ! Tperp Tperp = Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1 ! Parallel dynamic ! ddz derivative for Landau damping term Tpar = xnipp1j(ip) * ddz_nipj(ip+1,ij,iky,ikx,iz) & + xnipm1j(ip) * ddz_nipj(ip-1,ij,iky,ikx,iz) ! Mirror terms Tnipp1j = ynipp1j (ip,ij) * interp_nipj(ip+1,ij ,iky,ikx,iz) Tnipp1jm1 = ynipp1jm1(ip,ij) * interp_nipj(ip+1,ij-1,iky,ikx,iz) Tnipm1j = ynipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz) Tnipm1jm1 = ynipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz) ! Trapping terms Unipm1j = znipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz) Unipm1jp1 = znipm1jp1(ip,ij) * interp_nipj(ip-1,ij+1,iky,ikx,iz) Unipm1jm1 = znipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz) Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + Unipm1j + Unipm1jp1 + Unipm1jm1 !! Electrical potential term IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 Tphi = (xphij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) & + xphijp1_i(ip,ij)*kernel_i(ij+1,iky,ikx,iz,eo) & + xphijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*phikykxz ELSE Tphi = 0._dp ENDIF !! Vector potential term IF ( (p_int .LE. 3) .AND. (p_int .GE. 1) ) THEN ! Kronecker p1 or p3 Tpsi = (xpsij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) & + xpsijp1_i(ip,ij)*kernel_i(ij+1,iky,ikx,iz,eo) & + xpsijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*psikykxz ELSE Tpsi = 0._dp ENDIF !! Sum of all RHS terms moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = & ! Perpendicular magnetic gradient/curvature effects - -imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp & + -imagu*Ckxky(iky,ikx,iz,eo)*hatB(iz,eo) * Tperp & + ! Perpendicular pressure effects + -i_ky*beta*dpdx * Tperp& ! Parallel coupling (Landau damping) -gradz_coeff(iz,eo) * Tpar & ! Mirror term (parallel magnetic gradient) -gradzB(iz,eo)*gradz_coeff(iz,eo) * Tmir & ! Drives (density + temperature gradients) -i_ky * (Tphi - Tpsi) & ! Numerical hyperdiffusion (totally artificial, for stability purpose) -mu_x*diff_kx_coeff*kx**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) & -mu_y*diff_ky_coeff*ky**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) & ! Numerical parallel hyperdiffusion "mu_z*ddz**4" -mu_z*diff_dz_coeff*ddzND_nipj(ip,ij,iky,ikx,iz) & ! Collision term +TColl_i(ip,ij,iky,ikx,iz)& ! Nonlinear term with a (gxx*gxy - gxy**2)^1/2 factor -hatB_NL(iz,eo) * Sipj(ip,ij,iky,ikx,iz) IF(ip-4 .GT. 0) & ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33) moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = & moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) & + mu_p * moments_i(ip-4,ij,iky,ikx,iz,updatetlevel) ELSE moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp ENDIF END DO ploopi END DO jloopi END DO kyloopi END DO kxloopi END DO zloopi ! Execution time end CALL cpu_time(t1_rhs) tc_rhs = tc_rhs + (t1_rhs-t0_rhs) END SUBROUTINE moments_eq_rhs_i SUBROUTINE add_Maxwellian_background_terms ! This routine is meant to add the terms rising from the magnetic operator, ! i.e. (B x GradB) Grad, applied on the background Maxwellian distribution ! (x_a + spar^2)(b x GradB) GradFaM ! It gives birth to kx=ky=0 sources terms (averages) that hit moments 00, 20, ! 40, 01,02, 21 with background gradient dependences. USE prec_const USE time_integration, ONLY : updatetlevel USE model, ONLY: taue_qe, taui_qi, k_Ni, k_Ne, k_Ti, k_Te, KIN_E USE array, ONLY: moments_rhs_e, moments_rhs_i USE grid, ONLY: contains_kx0, contains_ky0, ikx_0, iky_0,& ips_e,ipe_e,ijs_e,ije_e,ips_i,ipe_i,ijs_i,ije_i,& jarray_e,parray_e,jarray_i,parray_i, zarray, izs,ize,& ip,ij IMPLICIT NONE real(dp), DIMENSION(izs:ize) :: sinz sinz(izs:ize) = SIN(zarray(izs:ize,0)) IF(contains_kx0 .AND. contains_ky0) THEN IF(KIN_E) THEN DO ip = ips_e,ipe_e DO ij = ijs_e,ije_e SELECT CASE(ij-1) CASE(0) ! j = 0 SELECT CASE (ip-1) CASE(0) ! Na00 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taue_qe * sinz(izs:ize) * (1.5_dp*k_Ne - 1.125_dp*k_Te) CASE(2) ! Na20 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*k_Ne - 2.75_dp*k_Te) CASE(4) ! Na40 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*k_Te END SELECT CASE(1) ! j = 1 SELECT CASE (ip-1) CASE(0) ! Na01 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& -taue_qe * sinz(izs:ize) * (k_Ne + 3.5_dp*k_Te) CASE(2) ! Na21 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& -taue_qe * sinz(izs:ize) * SQRT2*k_Te END SELECT CASE(2) ! j = 2 SELECT CASE (ip-1) CASE(0) ! Na02 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taue_qe * sinz(izs:ize) * 2._dp*k_Te END SELECT END SELECT ENDDO ENDDO ENDIF DO ip = ips_i,ipe_i DO ij = ijs_i,ije_i SELECT CASE(ij-1) CASE(0) ! j = 0 SELECT CASE (ip-1) CASE(0) ! Na00 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taui_qi * sinz(izs:ize) * (1.5_dp*k_Ni - 1.125_dp*k_Ti) CASE(2) ! Na20 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*k_Ni - 2.75_dp*k_Ti) CASE(4) ! Na40 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*k_Ti END SELECT CASE(1) ! j = 1 SELECT CASE (ip-1) CASE(0) ! Na01 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& -taui_qi * sinz(izs:ize) * (k_Ni + 3.5_dp*k_Ti) CASE(2) ! Na21 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& -taui_qi * sinz(izs:ize) * SQRT2*k_Ti END SELECT CASE(2) ! j = 2 SELECT CASE (ip-1) CASE(0) ! Na02 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taui_qi * sinz(izs:ize) * 2._dp*k_Ti END SELECT END SELECT ENDDO ENDDO ENDIF END SUBROUTINE END MODULE moments_eq_rhs diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90 index 8e9a2ce..6154717 100644 --- a/src/nonlinear_mod.F90 +++ b/src/nonlinear_mod.F90 @@ -1,463 +1,463 @@ MODULE nonlinear USE array, ONLY : dnjs, Sepj, Sipj, kernel_i, kernel_e,& moments_e_ZF, moments_i_ZF, phi_ZF USE initial_par, ONLY : ACT_ON_MODES USE basic USE fourier USE fields, ONLY : phi, psi, moments_e, moments_i USE grid USE model USE prec_const USE time_integration, ONLY : updatetlevel IMPLICIT NONE INCLUDE 'fftw3-mpi.f03' COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: F_cmpx, G_cmpx COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: Fx_cmpx, Gy_cmpx COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: Fy_cmpx, Gx_cmpx, F_conv_G INTEGER :: in, is, p_int, j_int INTEGER :: nmax, smax ! Upper bound of the sums REAL(dp):: kx, ky, kerneln, sqrt_p, sqrt_pp1 PUBLIC :: compute_Sapj, nonlinear_init CONTAINS SUBROUTINE nonlinear_init ALLOCATE( F_cmpx(ikys:ikye,ikxs:ikxe)) ALLOCATE( G_cmpx(ikys:ikye,ikxs:ikxe)) ALLOCATE(Fx_cmpx(ikys:ikye,ikxs:ikxe)) ALLOCATE(Gy_cmpx(ikys:ikye,ikxs:ikxe)) ALLOCATE(Fy_cmpx(ikys:ikye,ikxs:ikxe)) ALLOCATE(Gx_cmpx(ikys:ikye,ikxs:ikxe)) ALLOCATE(F_conv_G(ikys:ikye,ikxs:ikxe)) END SUBROUTINE nonlinear_init SUBROUTINE compute_Sapj ! This routine is meant to compute the non linear term for each specie and degree !! In real space Sapj ~ b*(grad(phi) x grad(g)) which in moments in fourier becomes !! Sapj = Sum_n (ikx Kn phi)#(iky Sum_s d_njs Naps) - (iky Kn phi)#(ikx Sum_s d_njs Naps) !! where # denotes the convolution. ! Execution time start CALL cpu_time(t0_Sapj) SELECT CASE(LINEARITY) CASE ('nonlinear') CALL compute_nonlinear CASE ('ZF_semilin') CALL compute_semi_linear_ZF CASE ('NZ_semilin') CALL compute_semi_linear_NZ CASE ('linear') Sepj = 0._dp; Sipj = 0._dp CASE DEFAULT IF(my_id.EQ.0) write(*,*) '/!\ Linearity not recognized /!\' stop END SELECT ! Execution time END CALL cpu_time(t1_Sapj) tc_Sapj = tc_Sapj + (t1_Sapj - t0_Sapj) END SUBROUTINE compute_Sapj SUBROUTINE compute_nonlinear IMPLICIT NONE !!!!!!!!!!!!!!!!!!!! ELECTRON non linear term computation (Sepj)!!!!!!!!!! IF(KIN_E) THEN zloope: DO iz = izs,ize ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments eo = MODULO(parray_e(ip),2) p_int = parray_e(ip) sqrt_p = SQRT(REAL(p_int,dp)); sqrt_pp1 = SQRT(REAL(p_int,dp)+1._dp); jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments j_int=jarray_e(ij) IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN !compute ! Set non linear sum truncation IF (NL_CLOS .EQ. -2) THEN nmax = Jmaxe ELSEIF (NL_CLOS .EQ. -1) THEN nmax = Jmaxe-j_int ELSE nmax = min(NL_CLOS,Jmaxe-j_int) ENDIF bracket_sum_r = 0._dp ! initialize sum over real nonlinear term nloope: DO in = 1,nmax+1 ! Loop over laguerre for the sum !-----------!! ELECTROSTATIC CONTRIBUTION {Sum_s dnjs Naps, Kernel phi} ! First convolution terms F_cmpx(ikys:ikye,ikxs:ikxe) = phi(ikys:ikye,ikxs:ikxe,iz) * kernel_e(in, ikys:ikye,ikxs:ikxe, iz, eo) ! Second convolution terms G_cmpx(ikys:ikye,ikxs:ikxe) = 0._dp ! initialization of the sum smax = MIN( (in-1)+(ij-1), Jmaxe ); DO is = 1, smax+1 ! sum truncation on number of moments G_cmpx(ikys:ikye,ikxs:ikxe) = G_cmpx(ikys:ikye,ikxs:ikxe) + & dnjs(in,ij,is) * moments_e(ip,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel) ENDDO !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\ CALL poisson_bracket_and_sum(F_cmpx,G_cmpx) !-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi} IF(EM) THEN ! First convolution terms F_cmpx(ikys:ikye,ikxs:ikxe) = -sqrt_tau_o_sigma_e * psi(ikys:ikye,ikxs:ikxe,iz) * kernel_e(in, ikys:ikye,ikxs:ikxe, iz, eo) ! Second convolution terms G_cmpx(ikys:ikye,ikxs:ikxe) = 0._dp ! initialization of the sum smax = MIN( (in-1)+(ij-1), Jmaxe ); DO is = 1, smax+1 ! sum truncation on number of moments G_cmpx(ikys:ikye,ikxs:ikxe) = G_cmpx(ikys:ikye,ikxs:ikxe) + & dnjs(in,ij,is) * (sqrt_pp1*moments_e(ip+1,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel)& +sqrt_p *moments_e(ip-1,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel)) ENDDO !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\ CALL poisson_bracket_and_sum(F_cmpx,G_cmpx) ENDIF ENDDO nloope !---------! Put back the real nonlinear product into k-space call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c) ! Retrieve convolution in input format DO iky = ikys, ikye Sepj(ip,ij,iky,ikxs:ikxe,iz) = bracket_sum_c(ikxs:ikxe,iky-local_nky_offset)*AA_x(ikxs:ikxe)*AA_y(iky) !Anti aliasing filter ENDDO ELSE Sepj(ip,ij,:,:,iz) = 0._dp ENDIF ENDDO jloope ENDDO ploope ENDDO zloope ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!! ION non linear term computation (Sipj)!!!!!!!!!! zloopi: DO iz = izs,ize ploopi: DO ip = ips_i,ipe_i ! Loop over Hermite moments p_int = parray_i(ip) eo = MODULO(parray_i(ip),2) jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments j_int=jarray_i(ij) - IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN !compute + IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN !compute for every moments except for closure 1 ! Set non linear sum truncation IF (NL_CLOS .EQ. -2) THEN nmax = Jmaxi ELSEIF (NL_CLOS .EQ. -1) THEN nmax = Jmaxi-j_int ELSE nmax = min(NL_CLOS,Jmaxi-j_int) ENDIF bracket_sum_r = 0._dp ! initialize sum over real nonlinear term nloopi: DO in = 1,nmax+1 ! Loop over laguerre for the sum -!-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi} +!-----------!! ELECTROSTATIC CONTRIBUTION ! First convolution terms F_cmpx(ikys:ikye,ikxs:ikxe) = phi(ikys:ikye,ikxs:ikxe,iz) * kernel_i(in, ikys:ikye,ikxs:ikxe, iz, eo) ! Second convolution terms G_cmpx(ikys:ikye,ikxs:ikxe) = 0._dp ! initialization of the sum smax = MIN( (in-1)+(ij-1), jmaxi ); DO is = 1, smax+1 ! sum truncation on number of moments G_cmpx(ikys:ikye,ikxs:ikxe) = G_cmpx(ikys:ikye,ikxs:ikxe) + & dnjs(in,ij,is) * moments_i(ip,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel) ENDDO !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\ CALL poisson_bracket_and_sum(F_cmpx,G_cmpx) !-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi} IF(EM) THEN ! First convolution terms F_cmpx(ikys:ikye,ikxs:ikxe) = -sqrt_tau_o_sigma_i * psi(ikys:ikye,ikxs:ikxe,iz) * kernel_i(in, ikys:ikye,ikxs:ikxe, iz, eo) ! Second convolution terms G_cmpx(ikys:ikye,ikxs:ikxe) = 0._dp ! initialization of the sum smax = MIN( (in-1)+(ij-1), Jmaxi ); DO is = 1, smax+1 ! sum truncation on number of moments G_cmpx(ikys:ikye,ikxs:ikxe) = G_cmpx(ikys:ikye,ikxs:ikxe) + & dnjs(in,ij,is) * (sqrt_pp1*moments_i(ip+1,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel)& +sqrt_p *moments_i(ip-1,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel)) ENDDO !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\ CALL poisson_bracket_and_sum(F_cmpx,G_cmpx) ENDIF ENDDO nloopi ! Put the real nonlinear product into k-space call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c) ! Retrieve convolution in input format DO iky = ikys, ikye Sipj(ip,ij,iky,ikxs:ikxe,iz) = bracket_sum_c(ikxs:ikxe,iky-local_nky_offset)*AA_x(ikxs:ikxe)*AA_y(iky) ENDDO ELSE Sipj(ip,ij,:,:,iz) = 0._dp ENDIF ENDDO jloopi ENDDO ploopi ENDDO zloopi !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END SUBROUTINE compute_nonlinear !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Semi linear computation : Only NZ-ZF convolutions are kept SUBROUTINE compute_semi_linear_ZF IMPLICIT NONE !!!!!!!!!!!!!!!!!!!! ELECTRON semi linear term computation (Sepj)!!!!!!!!!! IF(KIN_E) THEN zloope: DO iz = izs,ize ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments eo = MODULO(parray_e(ip),2) jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments j_int=jarray_e(ij) ! Set non linear sum truncation IF (NL_CLOS .EQ. -2) THEN nmax = Jmaxe ELSEIF (NL_CLOS .EQ. -1) THEN nmax = Jmaxe-(ij-1) ELSE nmax = NL_CLOS ENDIF bracket_sum_r = 0._dp ! initialize sum over real nonlinear term nloope: DO in = 1,nmax+1 ! Loop over laguerre for the sum ! Build the terms to convolve kxloope: DO ikx = ikxs,ikxe ! Loop over kx kyloope: DO iky = ikys,ikye ! Loop over ky kx = kxarray(ikx) ky = kyarray(iky) kerneln = kernel_e(in, ikx, iky, iz, eo) ! Zonal terms (=0 for all ky not 0) Fx_cmpx(iky,ikx) = 0._dp Gx_cmpx(iky,ikx) = 0._dp IF(iky .EQ. iky_0) THEN Fx_cmpx(iky,ikx) = imagu*kx* phi(iky,ikx,iz) * kerneln smax = MIN( (in-1)+(ij-1), jmaxe ); DO is = 1, smax+1 ! sum truncation on number of moments Gx_cmpx(iky,ikx) = Gx_cmpx(iky,ikx) + & dnjs(in,ij,is) * moments_e(ip,is,iky,ikx,iz,updatetlevel) ENDDO Gx_cmpx(iky,ikx) = imagu*kx*Gx_cmpx(iky,ikx) ENDIF ! NZ terms Fy_cmpx(iky,ikx) = imagu*ky* phi(iky,ikx,iz) * kerneln Gy_cmpx(iky,ikx) = 0._dp ! initialization of the sum smax = MIN( (in-1)+(ij-1), jmaxe ); DO is = 1, smax+1 ! sum truncation on number of moments Gy_cmpx(iky,ikx) = Gy_cmpx(iky,ikx) + & dnjs(in,ij,is) * moments_e(ip,is,iky,ikx,iz,updatetlevel) ENDDO Gy_cmpx(iky,ikx) = imagu*ky*Gy_cmpx(iky,ikx) ENDDO kyloope ENDDO kxloope ! First term df/dx x dg/dy CALL convolve_and_add(Fx_cmpx,Gy_cmpx) ! Second term -df/dy x dg/dx CALL convolve_and_add(-Fy_cmpx,Gx_cmpx) ENDDO nloope ! Put the real nonlinear product into k-space call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c) ! Retrieve convolution in input format DO ikx = ikxs, ikxe DO iky = ikys, ikye Sepj(ip,ij,iky,ikx,iz) = bracket_sum_c(ikx,iky-local_nky_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter ENDDO ENDDO ENDDO jloope ENDDO ploope ENDDO zloope ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!! ION non linear term computation (Sipj)!!!!!!!!!! zloopi: DO iz = izs,ize ploopi: DO ip = ips_i,ipe_i ! Loop over Hermite moments eo = MODULO(parray_i(ip),2) jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments j_int=jarray_i(ij) ! Set non linear sum truncation IF (NL_CLOS .EQ. -2) THEN nmax = Jmaxi ELSEIF (NL_CLOS .EQ. -1) THEN nmax = Jmaxi-(ij-1) ELSE nmax = NL_CLOS ENDIF bracket_sum_r = 0._dp ! initialize sum over real nonlinear term nloopi: DO in = 1,nmax+1 ! Loop over laguerre for the sum kxloopi: DO ikx = ikxs,ikxe ! Loop over kx kyloopi: DO iky = ikys,ikye ! Loop over ky ! Zonal terms (=0 for all ky not 0) Fx_cmpx(iky,ikx) = 0._dp Gx_cmpx(iky,ikx) = 0._dp IF(iky .EQ. iky_0) THEN Fx_cmpx(iky,ikx) = imagu*kx* phi(iky,ikx,iz) * kerneln smax = MIN( (in-1)+(ij-1), jmaxi ); DO is = 1, smax+1 ! sum truncation on number of moments Gx_cmpx(iky,ikx) = Gx_cmpx(iky,ikx) + & dnjs(in,ij,is) * moments_i(ip,is,iky,ikx,iz,updatetlevel) ENDDO Gx_cmpx(iky,ikx) = imagu*kx*Gx_cmpx(iky,ikx) ENDIF ! NZ terms Fy_cmpx(iky,ikx) = imagu*ky* phi(iky,ikx,iz) * kerneln Gy_cmpx(iky,ikx) = 0._dp ! initialization of the sum smax = MIN( (in-1)+(ij-1), jmaxi ); DO is = 1, smax+1 ! sum truncation on number of moments Gy_cmpx(iky,ikx) = Gy_cmpx(iky,ikx) + & dnjs(in,ij,is) * moments_i(ip,is,iky,ikx,iz,updatetlevel) ENDDO Gy_cmpx(iky,ikx) = imagu*ky*Gy_cmpx(iky,ikx) ENDDO kyloopi ENDDO kxloopi ! First term drphi x dzf CALL convolve_and_add(Fy_cmpx,Gx_cmpx) ! Second term -dzphi x drf CALL convolve_and_add(Fy_cmpx,Gx_cmpx) ENDDO nloopi ! Put the real nonlinear product into k-space call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c) ! Retrieve convolution in input format DO ikx = ikxs, ikxe DO iky = ikys, ikye Sipj(ip,ij,iky,ikx,iz) = bracket_sum_c(ikx,iky-local_nky_offset)*AA_x(ikx)*AA_y(iky) ENDDO ENDDO ENDDO jloopi ENDDO ploopi ENDDO zloopi END SUBROUTINE compute_semi_linear_ZF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Semi linear computation : Only kx=0*all convolutions are kept SUBROUTINE compute_semi_linear_NZ IMPLICIT NONE !!!!!!!!!!!!!!!!!!!! ELECTRON semi linear term computation (Sepj)!!!!!!!!!! IF(KIN_E) THEN zloope: DO iz = izs,ize ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments eo = MODULO(parray_e(ip),2) jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments j_int=jarray_e(ij) ! Set non linear sum truncation IF (NL_CLOS .EQ. -2) THEN nmax = Jmaxe ELSEIF (NL_CLOS .EQ. -1) THEN nmax = Jmaxe-(ij-1) ELSE nmax = NL_CLOS ENDIF bracket_sum_r = 0._dp ! initialize sum over real nonlinear term nloope: DO in = 1,nmax+1 ! Loop over laguerre for the sum ! Build the terms to convolve kxloope: DO ikx = ikxs,ikxe ! Loop over kx kyloope: DO iky = ikys,ikye ! Loop over ky kx = kxarray(ikx) ky = kyarray(iky) kerneln = kernel_e(in, ikx, iky, iz, eo) ! All terms Fx_cmpx(iky,ikx) = imagu*kx* phi(iky,ikx,iz) * kerneln smax = MIN( (in-1)+(ij-1), jmaxe ); DO is = 1, smax+1 ! sum truncation on number of moments Gx_cmpx(iky,ikx) = Gx_cmpx(iky,ikx) + & dnjs(in,ij,is) * moments_e(ip,is,iky,ikx,iz,updatetlevel) ENDDO Gx_cmpx(iky,ikx) = imagu*kx*Gx_cmpx(iky,ikx) ! Kx = 0 terms Fy_cmpx(iky,ikx) = 0._dp Gy_cmpx(iky,ikx) = 0._dp IF (ikx .EQ. ikx_0) THEN Fy_cmpx(iky,ikx) = imagu*ky* phi(iky,ikx,iz) * kerneln Gy_cmpx(iky,ikx) = 0._dp ! initialization of the sum smax = MIN( (in-1)+(ij-1), jmaxe ); DO is = 1, smax+1 ! sum truncation on number of moments Gy_cmpx(iky,ikx) = Gy_cmpx(iky,ikx) + & dnjs(in,ij,is) * moments_e(ip,is,iky,ikx,iz,updatetlevel) ENDDO Gy_cmpx(iky,ikx) = imagu*ky*Gy_cmpx(iky,ikx) ENDIF ENDDO kyloope ENDDO kxloope ! First term df/dx x dg/dy CALL convolve_and_add(Fx_cmpx,Gy_cmpx) ! Second term -df/dy x dg/dx CALL convolve_and_add(-Fy_cmpx,Gx_cmpx) ENDDO nloope ! Put the real nonlinear product into k-space call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c) ! Retrieve convolution in input format DO ikx = ikxs, ikxe DO iky = ikys, ikye Sepj(ip,ij,iky,ikx,iz) = bracket_sum_c(ikx,iky-local_nky_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter ENDDO ENDDO ENDDO jloope ENDDO ploope ENDDO zloope ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!! ION non linear term computation (Sipj)!!!!!!!!!! zloopi: DO iz = izs,ize ploopi: DO ip = ips_i,ipe_i ! Loop over Hermite moments eo = MODULO(parray_i(ip),2) jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments j_int=jarray_i(ij) ! Set non linear sum truncation IF (NL_CLOS .EQ. -2) THEN nmax = Jmaxi ELSEIF (NL_CLOS .EQ. -1) THEN nmax = Jmaxi-(ij-1) ELSE nmax = NL_CLOS ENDIF bracket_sum_r = 0._dp ! initialize sum over real nonlinear term nloopi: DO in = 1,nmax+1 ! Loop over laguerre for the sum kxloopi: DO ikx = ikxs,ikxe ! Loop over kx kyloopi: DO iky = ikys,ikye ! Loop over ky ! Zonal terms (=0 for all ky not 0) Fx_cmpx(iky,ikx) = imagu*kx* phi(iky,ikx,iz) * kerneln smax = MIN( (in-1)+(ij-1), jmaxi ); DO is = 1, smax+1 ! sum truncation on number of moments Gx_cmpx(iky,ikx) = Gx_cmpx(iky,ikx) + & dnjs(in,ij,is) * moments_i(ip,is,iky,ikx,iz,updatetlevel) ENDDO Gx_cmpx(iky,ikx) = imagu*kx*Gx_cmpx(iky,ikx) ! Kx = 0 terms Fy_cmpx(iky,ikx) = 0._dp Gy_cmpx(iky,ikx) = 0._dp IF (ikx .EQ. ikx_0) THEN Fy_cmpx(iky,ikx) = imagu*ky* phi(iky,ikx,iz) * kerneln Gy_cmpx(iky,ikx) = 0._dp ! initialization of the sum smax = MIN( (in-1)+(ij-1), jmaxi ); DO is = 1, smax+1 ! sum truncation on number of moments Gy_cmpx(iky,ikx) = Gy_cmpx(iky,ikx) + & dnjs(in,ij,is) * moments_i(ip,is,iky,ikx,iz,updatetlevel) ENDDO Gy_cmpx(iky,ikx) = imagu*ky*Gy_cmpx(iky,ikx) ENDIF ENDDO kyloopi ENDDO kxloopi ! First term drphi x dzf CALL convolve_and_add(Fy_cmpx,Gx_cmpx) ! Second term -dzphi x drf CALL convolve_and_add(Fy_cmpx,Gx_cmpx) ENDDO nloopi ! Put the real nonlinear product into k-space call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c) ! Retrieve convolution in input format DO ikx = ikxs, ikxe DO iky = ikys, ikye Sipj(ip,ij,iky,ikx,iz) = bracket_sum_c(ikx,iky-local_nky_offset)*AA_x(ikx)*AA_y(iky) ENDDO ENDDO ENDDO jloopi ENDDO ploopi ENDDO zloopi END SUBROUTINE compute_semi_linear_NZ END MODULE nonlinear diff --git a/src/processing_mod.F90 b/src/processing_mod.F90 index 54f99c6..7a4b278 100644 --- a/src/processing_mod.F90 +++ b/src/processing_mod.F90 @@ -1,800 +1,800 @@ MODULE processing USE basic USE prec_const USE grid implicit none REAL(dp), PUBLIC, PROTECTED :: pflux_ri, gflux_ri, pflux_re, gflux_re REAL(dp), PUBLIC, PROTECTED :: hflux_xi, hflux_xe PUBLIC :: compute_nadiab_moments_z_gradients_and_interp PUBLIC :: compute_density, compute_upar, compute_uperp PUBLIC :: compute_Tpar, compute_Tperp, compute_fluid_moments PUBLIC :: compute_radial_ion_transport, compute_radial_electron_transport PUBLIC :: compute_radial_ion_heatflux, compute_radial_electron_heatflux PUBLIC :: compute_Napjz_spectrum CONTAINS ! 1D diagnostic to compute the average radial particle transport _xyz SUBROUTINE compute_radial_ion_transport USE fields, ONLY : moments_i, phi, psi USE array, ONLY : kernel_i USE geometry, ONLY : Jacobian, iInt_Jacobian USE time_integration, ONLY : updatetlevel USE calculus, ONLY : simpson_rule_z USE model, ONLY : sqrt_tau_o_sigma_i, EM IMPLICIT NONE COMPLEX(dp) :: pflux_local, gflux_local, integral REAL(dp) :: ky_, buffer(1:2) INTEGER :: i_, world_rank, world_size, root COMPLEX(dp), DIMENSION(izgs:izge) :: integrant pflux_local = 0._dp ! particle flux gflux_local = 0._dp ! gyrocenter flux integrant = 0._dp ! auxiliary variable for z integration !!---------- Gyro center flux (drift kinetic) ------------ ! Electrostatic part IF(CONTAINS_ip0_i) THEN DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe integrant(izgs:izge) = integrant(izgs:izge) & +moments_i(ip0_i,ij0_i,iky,ikx,izgs:izge,updatetlevel)& *imagu*ky_*CONJG(phi(iky,ikx,izgs:izge)) ENDDO ENDDO ENDIF ! Electromagnetic part IF( EM .AND. CONTAINS_ip1_i ) THEN DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe integrant(izgs:izge) = integrant(izgs:izge)& -sqrt_tau_o_sigma_i*moments_i(ip1_i,ij0_i,iky,ikx,izgs:izge,updatetlevel)& *imagu*ky_*CONJG(psi(iky,ikx,izgs:izge)) ENDDO ENDDO ENDIF ! Integrate over z integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge) call simpson_rule_z(integrant,integral) ! Get process local gyrocenter flux gflux_local = integral*iInt_Jacobian ! integrant = 0._dp ! reset auxiliary variable !!---------- Particle flux (gyrokinetic) ------------ ! Electrostatic part IF(CONTAINS_ip0_i) THEN DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe DO ij = ijs_i, ije_i integrant(izgs:izge) = integrant(izgs:izge)& +moments_i(ip0_i,ij,iky,ikx,izgs:izge,updatetlevel)& *imagu*ky_*kernel_i(ij,iky,ikx,izgs:izge,0)*CONJG(phi(iky,ikx,izgs:izge)) ENDDO ENDDO ENDDO ENDIF ! Electromagnetic part IF( EM .AND. CONTAINS_ip1_i ) THEN DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe integrant = 0._dp ! auxiliary variable for z integration DO ij = ijs_i, ije_i integrant(izgs:izge) = integrant(izgs:izge)& -sqrt_tau_o_sigma_i*moments_i(ip1_i,ij,iky,ikx,izgs:izge,updatetlevel)& *imagu*ky_*kernel_i(ij,iky,ikx,izgs:izge,0)*CONJG(psi(iky,ikx,izgs:izge)) ENDDO ENDDO ENDDO ENDIF ! Integrate over z integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge) call simpson_rule_z(integrant,integral) ! Get process local particle flux pflux_local = integral*iInt_Jacobian !!!!---------- Sum over all processes ---------- buffer(1) = 2._dp*REAL(gflux_local) buffer(2) = 2._dp*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_ky .GT. 1) THEN !! Everyone sends its local_sum to root = 0 IF (rank_ky .NE. root) THEN CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr) ELSE ! Recieve from all the other processes DO i_ = 0,num_procs_ky-1 IF (i_ .NE. rank_ky) & CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr) gflux_ri = gflux_ri + buffer(1) pflux_ri = pflux_ri + buffer(2) ENDDO ENDIF ELSE gflux_ri = gflux_local pflux_ri = pflux_local ENDIF ! if(my_id .eq. 0) write(*,*) 'pflux_ri = ',pflux_ri END SUBROUTINE compute_radial_ion_transport ! 1D diagnostic to compute the average radial particle transport _xyz SUBROUTINE compute_radial_electron_transport USE fields, ONLY : moments_e, phi, psi USE array, ONLY : kernel_e USE geometry, ONLY : Jacobian, iInt_Jacobian USE time_integration, ONLY : updatetlevel USE calculus, ONLY : simpson_rule_z USE model, ONLY : sqrt_tau_o_sigma_e, EM IMPLICIT NONE COMPLEX(dp) :: pflux_local, gflux_local, integral REAL(dp) :: ky_, buffer(1:2) INTEGER :: i_, world_rank, world_size, root COMPLEX(dp), DIMENSION(izgs:izge) :: integrant pflux_local = 0._dp ! particle flux gflux_local = 0._dp ! gyrocenter flux integrant = 0._dp ! auxiliary variable for z integration !!---------- Gyro center flux (drift kinetic) ------------ ! Electrostatic part IF(CONTAINS_ip0_e) THEN DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe integrant(izgs:izge) = integrant(izgs:izge) & +moments_e(ip0_e,ij0_e,iky,ikx,izgs:izge,updatetlevel)& *imagu*ky_*CONJG(phi(iky,ikx,izgs:izge)) ENDDO ENDDO ENDIF ! Electromagnetic part IF( EM .AND. CONTAINS_ip1_e ) THEN DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe integrant(izgs:izge) = integrant(izgs:izge)& -sqrt_tau_o_sigma_e*moments_e(ip1_e,ij0_e,iky,ikx,izgs:izge,updatetlevel)& *imagu*ky_*CONJG(psi(iky,ikx,izgs:izge)) ENDDO ENDDO ENDIF ! Integrate over z integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge) call simpson_rule_z(integrant,integral) ! Get process local gyrocenter flux gflux_local = integral*iInt_Jacobian ! integrant = 0._dp ! reset auxiliary variable !!---------- Particle flux (gyrokinetic) ------------ ! Electrostatic part IF(CONTAINS_ip0_e) THEN DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe DO ij = ijs_e, ije_e integrant(izgs:izge) = integrant(izgs:izge)& +moments_e(ip0_e,ij,iky,ikx,izgs:izge,updatetlevel)& *imagu*ky_*kernel_e(ij,iky,ikx,izgs:izge,0)*CONJG(phi(iky,ikx,izgs:izge)) ENDDO ENDDO ENDDO ENDIF ! Electromagnetic part IF( EM .AND. CONTAINS_ip1_e ) THEN DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe integrant = 0._dp ! auxiliary variable for z integration DO ij = ijs_e, ije_e integrant(izgs:izge) = integrant(izgs:izge)& -sqrt_tau_o_sigma_e*moments_e(ip1_e,ij,iky,ikx,izgs:izge,updatetlevel)& *imagu*ky_*kernel_e(ij,iky,ikx,izgs:izge,0)*CONJG(psi(iky,ikx,izgs:izge)) ENDDO ENDDO ENDDO ENDIF ! Integrate over z integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge) call simpson_rule_z(integrant,integral) ! Get process local particle flux pflux_local = integral*iInt_Jacobian !!!!---------- Sum over all processes ---------- buffer(1) = 2._dp*REAL(gflux_local) buffer(2) = 2._dp*REAL(pflux_local) root = 0 !Gather manually among the rank_p=0 processes and perform the sum gflux_re = 0 pflux_re = 0 IF (num_procs_ky .GT. 1) THEN !! Everyone sends its local_sum to root = 0 IF (rank_ky .NE. root) THEN CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr) ELSE ! Recieve from all the other processes DO i_ = 0,num_procs_ky-1 IF (i_ .NE. rank_ky) & CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr) gflux_re = gflux_re + buffer(1) pflux_re = pflux_re + buffer(2) ENDDO ENDIF ELSE gflux_re = gflux_local pflux_re = pflux_local ENDIF END SUBROUTINE compute_radial_electron_transport ! 1D diagnostic to compute the average radial particle transport _xyz SUBROUTINE compute_radial_ion_heatflux USE fields, ONLY : moments_i, phi, psi USE array, ONLY : kernel_i, HF_phi_correction_operator USE geometry, ONLY : Jacobian, iInt_Jacobian USE time_integration, ONLY : updatetlevel USE model, ONLY : qi_taui, KIN_E, tau_i, EM USE calculus, ONLY : simpson_rule_z USE model, ONLY : sqrt_tau_o_sigma_i IMPLICIT NONE COMPLEX(dp) :: hflux_local, integral REAL(dp) :: ky_, buffer(1:2), n_dp INTEGER :: i_, world_rank, world_size, root, in COMPLEX(dp), DIMENSION(izgs:izge) :: integrant ! charge density q_a n_a hflux_local = 0._dp ! heat flux integrant = 0._dp ! z integration auxiliary variable !!----------------ELECTROSTATIC CONTRIBUTION--------------------------- IF(CONTAINS_ip0_i .AND. CONTAINS_ip2_i) THEN ! Loop to compute gamma_kx = sum_ky sum_j -i k_z Kernel_j Ni00 * phi DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe DO in = ijs_i, ije_i n_dp = jarray_i(in) integrant(izgs:izge) = integrant(izgs:izge) + tau_i*imagu*ky_*CONJG(phi(iky,ikx,izgs:izge))& *kernel_i(in,iky,ikx,izgs:izge,0)*(& 0.5_dp*SQRT2*moments_i(ip2_i,in ,iky,ikx,izgs:izge,updatetlevel)& +(2._dp*n_dp + 1.5_dp)*moments_i(ip0_i,in ,iky,ikx,izgs:izge,updatetlevel)& -(n_dp+1._dp)*moments_i(ip0_i,in+1,iky,ikx,izgs:izge,updatetlevel)& -n_dp*moments_i(ip0_i,in-1,iky,ikx,izgs:izge,updatetlevel)) ENDDO ENDDO ENDDO ENDIF IF(EM .AND. CONTAINS_ip1_i .AND. CONTAINS_ip3_i) THEN !!----------------ELECTROMAGNETIC CONTRIBUTION-------------------- DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe DO in = ijs_i, ije_i n_dp = jarray_i(in) integrant(izgs:izge) = integrant(izgs:izge) & +tau_i*sqrt_tau_o_sigma_i*imagu*ky_*CONJG(psi(iky,ikx,izgs:izge))& *kernel_i(in,iky,ikx,izgs:izge,0)*(& 0.5_dp*SQRT2*SQRT3*moments_i(ip3_i,in ,iky,ikx,izgs:izge,updatetlevel)& +1.5_dp*moments_i(ip1_i,in ,iky,ikx,izgs:izge,updatetlevel)& +(2._dp*n_dp+1._dp)*moments_i(ip1_i,in ,iky,ikx,izgs:izge,updatetlevel)& -(n_dp+1._dp)*moments_i(ip1_i,in+1,iky,ikx,izgs:izge,updatetlevel)& -n_dp*moments_i(ip1_i,in-1,iky,ikx,izgs:izge,updatetlevel)) ENDDO ENDDO ENDDO ENDIF ! Add polarisation contribution - integrant(izs:ize) = integrant(izs:ize) + tau_i*imagu*ky_& - *CONJG(phi(iky,ikx,izs:ize))*phi(iky,ikx,izs:ize) * HF_phi_correction_operator(iky,ikx,izs:ize) + ! integrant(izgs:izge) = integrant(izgs:izge) + tau_i*imagu*ky_& + ! *CONJG(phi(iky,ikx,izgs:izge))*phi(iky,ikx,izgs:izge) * HF_phi_correction_operator(iky,ikx,izgs:izge) ! Integrate over z integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge) call simpson_rule_z(integrant,integral) hflux_local = hflux_local + integral*iInt_Jacobian ! Double it for taking into account the other half plane buffer(2) = 2._dp*REAL(hflux_local) root = 0 !Gather manually among the rank_p=0 processes and perform the sum hflux_xi = 0 IF (num_procs_ky .GT. 1) THEN !! Everyone sends its local_sum to root = 0 IF (rank_ky .NE. root) THEN CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr) ELSE ! Recieve from all the other processes DO i_ = 0,num_procs_ky-1 IF (i_ .NE. rank_ky) & CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr) hflux_xi = hflux_xi + buffer(2) ENDDO ENDIF ELSE hflux_xi = hflux_local ENDIF END SUBROUTINE compute_radial_ion_heatflux ! 1D diagnostic to compute the average radial particle transport _xyz SUBROUTINE compute_radial_electron_heatflux USE fields, ONLY : moments_e, phi, psi USE array, ONLY : kernel_e, HF_phi_correction_operator USE geometry, ONLY : Jacobian, iInt_Jacobian USE time_integration, ONLY : updatetlevel USE model, ONLY : qi_taui, KIN_E, tau_e, EM USE calculus, ONLY : simpson_rule_z USE model, ONLY : sqrt_tau_o_sigma_e IMPLICIT NONE COMPLEX(dp) :: hflux_local, integral REAL(dp) :: ky_, buffer(1:2), n_dp INTEGER :: i_, world_rank, world_size, root, in COMPLEX(dp), DIMENSION(izgs:izge) :: integrant ! charge density q_a n_a hflux_local = 0._dp ! heat flux integrant = 0._dp ! z integration auxiliary variable !!----------------ELECTROSTATIC CONTRIBUTION--------------------------- IF(CONTAINS_ip0_e .AND. CONTAINS_ip2_e) THEN ! Loop to compute gamma_kx = sum_ky sum_j -i k_z Kernel_j Ni00 * phi DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe DO in = ijs_e, ije_e n_dp = jarray_e(in) integrant(izgs:izge) = integrant(izgs:izge) + tau_e*imagu*ky_*CONJG(phi(iky,ikx,izgs:izge))& *kernel_e(in,iky,ikx,izgs:izge,0)*(& 0.5_dp*SQRT2*moments_e(ip2_e,in ,iky,ikx,izgs:izge,updatetlevel)& +(2._dp*n_dp + 1.5_dp)*moments_e(ip0_e,in ,iky,ikx,izgs:izge,updatetlevel)& -(n_dp+1._dp)*moments_e(ip0_e,in+1,iky,ikx,izgs:izge,updatetlevel)& -n_dp*moments_e(ip0_e,in-1,iky,ikx,izgs:izge,updatetlevel)) ENDDO ENDDO ENDDO ENDIF IF(EM .AND. CONTAINS_ip1_e .AND. CONTAINS_ip3_e) THEN !!----------------ELECTROMAGNETIC CONTRIBUTION-------------------- DO iky = ikys,ikye ky_ = kyarray(iky) DO ikx = ikxs,ikxe DO in = ijs_e, ije_e n_dp = jarray_e(in) integrant(izgs:izge) = integrant(izgs:izge) & +tau_e*sqrt_tau_o_sigma_e*imagu*ky_*CONJG(psi(iky,ikx,izgs:izge))& *kernel_e(in,iky,ikx,izgs:izge,0)*(& 0.5_dp*SQRT2*SQRT3*moments_e(ip3_e,in ,iky,ikx,izgs:izge,updatetlevel)& +1.5_dp*CONJG(moments_e(ip1_e,in ,iky,ikx,izgs:izge,updatetlevel))& !????? +(2._dp*n_dp+1._dp)*moments_e(ip1_e,in ,iky,ikx,izgs:izge,updatetlevel)& -(n_dp+1._dp)*moments_e(ip1_e,in+1,iky,ikx,izgs:izge,updatetlevel)& -n_dp*moments_e(ip1_e,in-1,iky,ikx,izgs:izge,updatetlevel)) ENDDO ENDDO ENDDO ENDIF ! Add polarisation contribution integrant(izs:ize) = integrant(izs:ize) + tau_e*imagu*ky_& *CONJG(phi(iky,ikx,izs:ize))*phi(iky,ikx,izs:ize) * HF_phi_correction_operator(iky,ikx,izs:ize) ! Integrate over z integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge) call simpson_rule_z(integrant,integral) hflux_local = hflux_local + integral*iInt_Jacobian ! Double it for taking into account the other half plane buffer(2) = 2._dp*REAL(hflux_local) root = 0 !Gather manually among the rank_p=0 processes and perform the sum hflux_xi = 0 IF (num_procs_ky .GT. 1) THEN !! Everyone sends its local_sum to root = 0 IF (rank_ky .NE. root) THEN CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr) ELSE ! Recieve from all the other processes DO i_ = 0,num_procs_ky-1 IF (i_ .NE. rank_ky) & CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr) hflux_xi = hflux_xi + buffer(2) ENDDO ENDIF ELSE hflux_xi = hflux_local ENDIF END SUBROUTINE compute_radial_electron_heatflux SUBROUTINE compute_nadiab_moments_z_gradients_and_interp ! evaluate the non-adiabatique ion moments ! ! n_{pi} = N^{pj} + kernel(j) /tau_i phi delta_p0 ! USE fields, ONLY : moments_i, moments_e, phi, psi USE array, ONLY : kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i, & ddz_nepj, ddzND_nepj, interp_nepj,& ddz_nipj, ddzND_nipj, interp_nipj USE time_integration, ONLY : updatetlevel USE model, ONLY : qe_taue, qi_taui,q_o_sqrt_tau_sigma_e, q_o_sqrt_tau_sigma_i, & KIN_E, CLOS, beta USE calculus, ONLY : grad_z, grad_z2, grad_z4, interp_z IMPLICIT NONE INTEGER :: eo, p_int, j_int CALL cpu_time(t0_process) ! Electron non adiab moments IF(KIN_E) THEN DO ip=ipgs_e,ipge_e IF(parray_e(ip) .EQ. 0) THEN DO ij=ijgs_e,ijge_e nadiab_moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge) = moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) & + qe_taue*kernel_e(ij,ikys:ikye,ikxs:ikxe,izgs:izge,0)*phi(ikys:ikye,ikxs:ikxe,izgs:izge) ENDDO ELSEIF( (parray_e(ip) .EQ. 1) .AND. (beta .GT. 0) ) THEN DO ij=ijgs_e,ijge_e nadiab_moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge) = moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) & - q_o_sqrt_tau_sigma_e*kernel_e(ij,ikys:ikye,ikxs:ikxe,izgs:izge,0)*psi(ikys:ikye,ikxs:ikxe,izgs:izge) ENDDO ELSE DO ij=ijgs_e,ijge_e nadiab_moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge) = moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) ENDDO ENDIF ENDDO ENDIF ! Ions non adiab moments DO ip=ipgs_i,ipge_i IF(parray_i(ip) .EQ. 0) THEN DO ij=ijgs_i,ijge_i nadiab_moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge) = moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) & + qi_taui*kernel_i(ij,ikys:ikye,ikxs:ikxe,izgs:izge,0)*phi(ikys:ikye,ikxs:ikxe,izgs:izge) ENDDO ELSEIF( (parray_i(ip) .EQ. 1) .AND. (beta .GT. 0) ) THEN DO ij=ijgs_i,ijge_i nadiab_moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge) = moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) & - q_o_sqrt_tau_sigma_i*kernel_i(ij,ikys:ikye,ikxs:ikxe,izgs:izge,0)*psi(ikys:ikye,ikxs:ikxe,izgs:izge) ENDDO ELSE DO ij=ijgs_i,ijge_i nadiab_moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge) = moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) ENDDO ENDIF ENDDO !! Ensure to kill the moments too high if closue option is set to 1 IF(CLOS .EQ. 1) THEN IF(KIN_E) THEN DO ip=ipgs_e,ipge_e p_int = parray_e(ip) DO ij=ijgs_e,ijge_e j_int = jarray_e(ij) IF(p_int+2*j_int .GT. dmaxe) & nadiab_moments_e(ip,ij,:,:,:) = 0._dp ENDDO ENDDO ENDIF DO ip=ipgs_i,ipge_i p_int = parray_i(ip) DO ij=ijgs_i,ijge_i j_int = jarray_i(ij) IF(p_int+2*j_int .GT. dmaxi) & nadiab_moments_i(ip,ij,:,:,:) = 0._dp ENDDO ENDDO ENDIF !------------- INTERP AND GRADIENTS ALONG Z ---------------------------------- IF (KIN_E) THEN DO ikx = ikxs,ikxe DO iky = ikys,ikye DO ij = ijgs_e,ijge_e DO ip = ipgs_e,ipge_e p_int = parray_e(ip) eo = MODULO(p_int,2) ! Indicates if we are on even or odd z grid ! Compute z derivatives CALL grad_z(eo,nadiab_moments_e(ip,ij,iky,ikx,izgs:izge), ddz_nepj(ip,ij,iky,ikx,izs:ize)) ! Parallel hyperdiffusion CALL grad_z4(moments_e(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddzND_nepj(ip,ij,iky,ikx,izs:ize)) ! CALL grad_z2(nadiab_moments_e(ip,ij,iky,ikx,izgs:izge),ddzND_nepj(ip,ij,iky,ikx,izs:ize)) ! Compute even odd grids interpolation CALL interp_z(eo,nadiab_moments_e(ip,ij,iky,ikx,izgs:izge), interp_nepj(ip,ij,iky,ikx,izs:ize)) ENDDO ENDDO ENDDO ENDDO ENDIF DO ikx = ikxs,ikxe DO iky = ikys,ikye DO ij = ijgs_i,ijge_i DO ip = ipgs_i,ipge_i p_int = parray_i(ip) eo = MODULO(p_int,2) ! Indicates if we are on even or odd z grid ! Compute z first derivative CALL grad_z(eo,nadiab_moments_i(ip,ij,iky,ikx,izgs:izge), ddz_nipj(ip,ij,iky,ikx,izs:ize)) ! Parallel numerical diffusion CALL grad_z4(moments_i(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddzND_nipj(ip,ij,iky,ikx,izs:ize)) ! CALL grad_z2(nadiab_moments_i(ip,ij,iky,ikx,izgs:izge),ddzND_nipj(ip,ij,iky,ikx,izs:ize)) ! Compute even odd grids interpolation CALL interp_z(eo,nadiab_moments_i(ip,ij,iky,ikx,izgs:izge), interp_nipj(ip,ij,iky,ikx,izs:ize)) ENDDO ENDDO ENDDO ENDDO ! Execution time end CALL cpu_time(t1_process) tc_process = tc_process + (t1_process - t0_process) END SUBROUTINE compute_nadiab_moments_z_gradients_and_interp SUBROUTINE compute_Napjz_spectrum USE fields, ONLY : moments_e, moments_i USE model, ONLY : KIN_E USE array, ONLY : Nipjz, Nepjz USE time_integration, ONLY : updatetlevel IMPLICIT NONE REAL(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize) :: local_sum_e,global_sum_e, buffer_e REAL(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize) :: local_sum_i,global_sum_i, buffer_i INTEGER :: i_, world_rank, world_size, root, count root = 0 ! Electron moments spectrum IF (KIN_E) THEN ! build local sum local_sum_e = 0._dp DO ikx = ikxs,ikxe DO iky = ikys,ikye local_sum_e(:,:,:) = local_sum_e(:,:,:) + & moments_e(:,:,iky,ikx,:,updatetlevel) * CONJG(moments_e(:,:,iky,ikx,:,updatetlevel)) ENDDO ENDDO ! sum reduction buffer_e = local_sum_e global_sum_e = 0._dp count = (ipe_e-ips_e+1)*(ije_e-ijs_e+1)*(ize-izs+1) IF (num_procs_ky .GT. 1) THEN !! Everyone sends its local_sum to root = 0 IF (rank_ky .NE. root) THEN CALL MPI_SEND(buffer_e, count , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr) ELSE ! Recieve from all the other processes DO i_ = 0,num_procs_ky-1 IF (i_ .NE. rank_ky) & CALL MPI_RECV(buffer_e, count , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr) global_sum_e = global_sum_e + buffer_e ENDDO ENDIF ELSE global_sum_e = local_sum_e ENDIF Nepjz = global_sum_e ENDIF ! Ion moment spectrum ! build local sum local_sum_i = 0._dp DO ikx = ikxs,ikxe DO iky = ikys,ikye local_sum_i(:,:,:) = local_sum_i(:,:,:) + & moments_i(:,:,iky,ikx,:,updatetlevel) * CONJG(moments_i(:,:,iky,ikx,:,updatetlevel)) ENDDO ENDDO ! sum reduction buffer_i = local_sum_i global_sum_i = 0._dp count = (ipe_i-ips_i+1)*(ije_i-ijs_i+1)*(ize-izs+1) IF (num_procs_ky .GT. 1) THEN !! Everyone sends its local_sum to root = 0 IF (rank_ky .NE. root) THEN CALL MPI_SEND(buffer_i, count , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr) ELSE ! Recieve from all the other processes DO i_ = 0,num_procs_ky-1 IF (i_ .NE. rank_ky) & CALL MPI_RECV(buffer_i, count , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr) global_sum_i = global_sum_i + buffer_i ENDDO ENDIF ELSE global_sum_i = local_sum_i ENDIF Nipjz = global_sum_i END SUBROUTINE compute_Napjz_spectrum !_____________________________________________________________________________! !!!!! FLUID MOMENTS COMPUTATIONS !!!!! ! Compute the 2D particle density for electron and ions (sum over Laguerre) SUBROUTINE compute_density USE array, ONLY : dens_e, dens_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i USE model, ONLY : KIN_E IMPLICIT NONE COMPLEX(dp) :: dens IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i ) THEN ! Loop to compute dens_i = sum_j kernel_j Ni0j at each k DO iz = izs,ize DO iky = ikys,ikye DO ikx = ikxs,ikxe IF(KIN_E) THEN ! electron density dens = 0._dp DO ij = ijs_e, ije_e dens = dens + kernel_e(ij,iky,ikx,iz,0) * nadiab_moments_e(ip0_e,ij,iky,ikx,iz) ENDDO dens_e(iky,ikx,iz) = dens ENDIF ! ion density dens = 0._dp DO ij = ijs_i, ije_i dens = dens + kernel_i(ij,iky,ikx,iz,0) * nadiab_moments_i(ip0_e,ij,iky,ikx,iz) ENDDO dens_i(iky,ikx,iz) = dens ENDDO ENDDO ENDDO ENDIF END SUBROUTINE compute_density ! Compute the 2D particle fluid perp velocity for electron and ions (sum over Laguerre) SUBROUTINE compute_uperp USE array, ONLY : uper_e, uper_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i USE model, ONLY : KIN_E IMPLICIT NONE COMPLEX(dp) :: uperp IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i ) THEN DO iz = izs,ize DO iky = ikys,ikye DO ikx = ikxs,ikxe IF(KIN_E) THEN ! electron uperp = 0._dp DO ij = ijs_e, ije_e uperp = uperp + kernel_e(ij,iky,ikx,iz,0) *& 0.5_dp*(nadiab_moments_e(ip0_e,ij,iky,ikx,iz) - nadiab_moments_e(ip0_e,ij-1,iky,ikx,iz)) ENDDO uper_e(iky,ikx,iz) = uperp ENDIF ! ion uperp = 0._dp DO ij = ijs_i, ije_i uperp = uperp + kernel_i(ij,iky,ikx,iz,0) *& 0.5_dp*(nadiab_moments_i(ip0_i,ij,iky,ikx,iz) - nadiab_moments_i(ip0_i,ij-1,iky,ikx,iz)) ENDDO uper_i(iky,ikx,iz) = uperp ENDDO ENDDO ENDDO ENDIF END SUBROUTINE compute_uperp ! Compute the 2D particle fluid par velocity for electron and ions (sum over Laguerre) SUBROUTINE compute_upar USE array, ONLY : upar_e, upar_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i USE model, ONLY : KIN_E IMPLICIT NONE COMPLEX(dp) :: upar IF ( CONTAINS_ip1_e .AND. CONTAINS_ip1_i ) THEN DO iz = izs,ize DO iky = ikys,ikye DO ikx = ikxs,ikxe IF(KIN_E) THEN ! electron upar = 0._dp DO ij = ijs_e, ije_e upar = upar + kernel_e(ij,iky,ikx,iz,1)*nadiab_moments_e(ip1_e,ij,iky,ikx,iz) ENDDO upar_e(iky,ikx,iz) = upar ENDIF ! ion upar = 0._dp DO ij = ijs_i, ije_i upar = upar + kernel_i(ij,iky,ikx,iz,1)*nadiab_moments_i(ip1_i,ij,iky,ikx,iz) ENDDO upar_i(iky,ikx,iz) = upar ENDDO ENDDO ENDDO ELSE IF(KIN_E)& upar_e = 0 upar_i = 0 ENDIF END SUBROUTINE compute_upar ! Compute the 2D particle temperature for electron and ions (sum over Laguerre) SUBROUTINE compute_tperp USE array, ONLY : Tper_e, Tper_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i USE model, ONLY : KIN_E IMPLICIT NONE REAL(dp) :: j_dp COMPLEX(dp) :: Tperp IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i .AND. & CONTAINS_ip2_e .AND. CONTAINS_ip2_i ) THEN ! Loop to compute T = 1/3*(Tpar + 2Tperp) DO iz = izs,ize DO iky = ikys,ikye DO ikx = ikxs,ikxe ! electron temperature IF(KIN_E) THEN Tperp = 0._dp DO ij = ijs_e, ije_e j_dp = REAL(ij-1,dp) Tperp = Tperp + kernel_e(ij,iky,ikx,iz,0)*& ((2_dp*j_dp+1)*nadiab_moments_e(ip0_e,ij ,iky,ikx,iz)& -j_dp *nadiab_moments_e(ip0_e,ij-1,iky,ikx,iz)& -j_dp+1 *nadiab_moments_e(ip0_e,ij+1,iky,ikx,iz)) ENDDO Tper_e(iky,ikx,iz) = Tperp ENDIF ! ion temperature Tperp = 0._dp DO ij = ijs_i, ije_i j_dp = REAL(ij-1,dp) Tperp = Tperp + kernel_i(ij,iky,ikx,iz,0)*& ((2_dp*j_dp+1)*nadiab_moments_i(ip0_i,ij ,iky,ikx,iz)& -j_dp *nadiab_moments_i(ip0_i,ij-1,iky,ikx,iz)& -j_dp+1 *nadiab_moments_i(ip0_i,ij+1,iky,ikx,iz)) ENDDO Tper_i(iky,ikx,iz) = Tperp ENDDO ENDDO ENDDO ENDIF END SUBROUTINE compute_Tperp ! Compute the 2D particle temperature for electron and ions (sum over Laguerre) SUBROUTINE compute_Tpar USE array, ONLY : Tpar_e, Tpar_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i USE model, ONLY : KIN_E IMPLICIT NONE REAL(dp) :: j_dp COMPLEX(dp) :: tpar IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i .AND. & CONTAINS_ip2_e .AND. CONTAINS_ip2_i ) THEN ! Loop to compute T = 1/3*(Tpar + 2Tperp) DO iz = izs,ize DO iky = ikys,ikye DO ikx = ikxs,ikxe ! electron temperature IF(KIN_E) THEN Tpar = 0._dp DO ij = ijs_e, ije_e j_dp = REAL(ij-1,dp) Tpar = Tpar + kernel_e(ij,iky,ikx,iz,0)*& (SQRT2 * nadiab_moments_e(ip2_e,ij,iky,ikx,iz) & + nadiab_moments_e(ip0_e,ij,iky,ikx,iz)) ENDDO Tpar_e(iky,ikx,iz) = Tpar ENDIF ! ion temperature Tpar = 0._dp DO ij = ijs_i, ije_i j_dp = REAL(ij-1,dp) Tpar = Tpar + kernel_i(ij,iky,ikx,iz,0)*& (SQRT2 * nadiab_moments_i(ip2_i,ij,iky,ikx,iz) & + nadiab_moments_i(ip0_i,ij,iky,ikx,iz)) ENDDO Tpar_i(iky,ikx,iz) = Tpar ENDDO ENDDO ENDDO ENDIF END SUBROUTINE compute_Tpar ! Compute the 2D particle fluid moments for electron and ions (sum over Laguerre) SUBROUTINE compute_fluid_moments USE array, ONLY : dens_e, Tpar_e, Tper_e, dens_i, Tpar_i, Tper_i, temp_e, temp_i USE model, ONLY : KIN_E IMPLICIT NONE CALL compute_density CALL compute_upar CALL compute_uperp CALL compute_Tpar CALL compute_Tperp ! Temperature temp_e = (Tpar_e - 2._dp * Tper_e)/3._dp - dens_e temp_i = (Tpar_i - 2._dp * Tper_i)/3._dp - dens_i END SUBROUTINE compute_fluid_moments END MODULE processing diff --git a/testcases/miller_example/fort_00.90 b/testcases/miller_example/bug/fort_00.90 similarity index 93% copy from testcases/miller_example/fort_00.90 copy to testcases/miller_example/bug/fort_00.90 index 6958e7a..728a345 100644 Binary files a/testcases/miller_example/fort_00.90 and b/testcases/miller_example/bug/fort_00.90 differ diff --git a/testcases/miller_example/fort_00.90 b/testcases/miller_example/bug/fort_01.90 similarity index 90% copy from testcases/miller_example/fort_00.90 copy to testcases/miller_example/bug/fort_01.90 index 6958e7a..7deb16c 100644 Binary files a/testcases/miller_example/fort_00.90 and b/testcases/miller_example/bug/fort_01.90 differ diff --git a/testcases/miller_example/fort_00.90 b/testcases/miller_example/fort_00.90 index 6958e7a..728a345 100644 Binary files a/testcases/miller_example/fort_00.90 and b/testcases/miller_example/fort_00.90 differ diff --git a/testcases/miller_example/fort_00.90 b/testcases/miller_example_w_salpha/fort_00.90 similarity index 91% copy from testcases/miller_example/fort_00.90 copy to testcases/miller_example_w_salpha/fort_00.90 index 6958e7a..2701de9 100644 Binary files a/testcases/miller_example/fort_00.90 and b/testcases/miller_example_w_salpha/fort_00.90 differ diff --git a/testcases/miller_example/fort_00.90 b/testcases/miller_example_w_salpha/fort_01.90 similarity index 90% copy from testcases/miller_example/fort_00.90 copy to testcases/miller_example_w_salpha/fort_01.90 index 6958e7a..7deb16c 100644 Binary files a/testcases/miller_example/fort_00.90 and b/testcases/miller_example_w_salpha/fort_01.90 differ diff --git a/wk/Zpinch_coll_scan_kN_1.7.m b/wk/Zpinch_coll_scan_kN_1.7.m index e223e0e..ad00ba1 100644 --- a/wk/Zpinch_coll_scan_kN_1.7.m +++ b/wk/Zpinch_coll_scan_kN_1.7.m @@ -1,68 +1,83 @@ -%% if 0 +%% figure Kn = 1.7; % SUGAMA DK 4,2 % nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00]; % Gavg_a = 1e-2*[1.00 1.71 2.18 3.11 4.11 5.20 6.08]; % Gstd_a = 1e-2*[1.78 2.67 2.82 3.08 2.33 1.35 1.43]; % errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama DK (4,2)'); hold on -% SUGAMA GK 4,2 +% SUGAMA GK 4,2 200x32 nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.0]; -Gavg_a = [2.54e-2 4.66e-2 6.96e-2 8.98e-2 1.06e-1 1.24e-1 1.43e-1 1.52e-1 1.69e-1 1.09e-1]; +Gavg_a = [2.54e-2 4.66e-2 6.96e-2 8.98e-2 1.06e-1 1.24e-1 1.43e-1 1.52e-1 1.69e-1 1.79e-1]; Gstd_a = [3.04e-2 1.42e-2 1.56e-2 1.23e-2 1.20e-2 1.57e-2 1.63e-2 2.06e-2 2.14e-02 1.78e-2]; -errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama GK (4,2)'); hold on +errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama GK (4,2) 200x32'); hold on + +% SUGAMA GK 4,2 256x64 +nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.0]; +Gavg_a = [1.99e-2 3.03e-2 4.44e-2 5.87e-2 7.87e-2 1.07e-1 1.22e-1 1.38e-1 1.63e-1 1.75e-1]; +Gstd_a = [1.55e-2 2.22e-2 1.12e-2 1.54e-2 2.20e-2 2.76e-2 2.44e-2 3.20e-2 3.12e-2 3.13e-2]; +errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama GK (4,2) 256x64'); hold on -% FCGK 4,2 + +% SUGAMA GK 6,3 200x32 +nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.0]; +Gavg_a = [2.35e-2 3.57e-2 5.09e-2 6.97e-2 8.65e-2 1.01e-1 1.16e-1 1.33e-1 1.49e-1 1.62e-1]; +Gstd_a = [3.76e-2 3.91e-2 2.96e-2 1.57e-2 1.73e-2 1.94e-2 2.21e-2 2.01e-2 1.93e-2 2.24e-2]; + +errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama GK (6,3) 200x32'); hold on + + +% FCGK 4,2 200x32 nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 10.0]; Gavg_a = [8.57e-2 1.45e-1 2.25e-1 2.87e-1 3.48e-1 4.06e-1 4.51e-1 3.65e-1*Kn]; Gstd_a = [2.07e-2 2.61e-2 2.40e-2 3.46e-2 4.30e-2 5.00e-2 5.11e-2 0]; -errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Coulomb (4,2)'); hold on - -% LDGK ii 6,3 -nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00]; -Gavg_a = [3.86e-2 1.82e-2 3.08e-2 5.24e-2 7.08e-2 8.26e-2 5.78e-2 7.16e-2 7.96e-2]; -Gstd_a = [3.52e-2 1.87e-2 2.86e-2 2.79e-2 1.72e-2 2.40e-2 2.46e-2 1.01e-2 1.21e-2]; +errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Coulomb (4,2) 200x32'); hold on -errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Landau ii (6,3)'); hold on +% % LDGK ii 6,3 200x32 +% nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00]; +% Gavg_a = [3.86e-2 1.82e-2 3.08e-2 5.24e-2 7.08e-2 8.26e-2 5.78e-2 7.16e-2 7.96e-2]; +% Gstd_a = [3.52e-2 1.87e-2 2.86e-2 2.79e-2 1.72e-2 2.40e-2 2.46e-2 1.01e-2 1.21e-2]; +% +% errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Landau ii (6,3) 200x32'); hold on % Collisionless plot([0 1], 0.02343*[1 1],'--k','DisplayName','$\nu=0$'); % xlim([0 0.1]); legend('show'); xlabel('$\nu R/c_s$'); ylabel('$\Gamma_x^\infty/\kappa_N$'); - +set(gca,'Yscale','log'); end if 0 %% figure nu = 0.1; % FCGK 4,2 kn_a = [1.60 1.80 2.00 2.20 2.40]; Gavg_a = [1.11e-1 6.86e-1 3.44e-0 1.12e+1 2.87e+1]; Gstd_a = [7.98e-3 1.10e-1 4.03e-1 2.03e+0 7.36e+0]; errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Coulomb (4,2)'); hold on % % Collisionless % plot([0 1], 0.02343*[1 1],'--k','DisplayName','$\nu=0$'); % xlim([1.6 2.5]); legend('show'); xlabel('$\nu R/c_s$'); ylabel('$\Gamma_x^\infty/\kappa_N$'); end \ No newline at end of file diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m index f45fc16..de09d16 100644 --- a/wk/analysis_gene.m +++ b/wk/analysis_gene.m @@ -1,158 +1,159 @@ helazdir = '/home/ahoffman/HeLaZ/'; addpath(genpath([helazdir,'matlab'])) % ... add addpath(genpath([helazdir,'matlab/plot'])) % ... add addpath(genpath([helazdir,'matlab/compute'])) % ... add addpath(genpath([helazdir,'matlab/load'])) % ... add % folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/'; % folder = '/misc/gene_results/shearless_cyclone/miller_output_0.8/'; % folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/'; % folder = '/misc/gene_results/shearless_cyclone/rm_corrections_HF/'; % folder = '/misc/gene_results/shearless_cyclone/linear_s_alpha_CBC_100/'; % folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/'; % folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_1.0/'; % folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_0.8/'; % folder = '/misc/gene_results/HP_fig_2a_mu_1e-2/'; % folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/'; % folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/'; % folder = '/misc/gene_results/LD_zpinch_1.6/'; % folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/'; % folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/'; % folder = '/misc/gene_results/Z-pinch/ZP_HP_kn_1.6_HRES/'; % folder = '/misc/gene_results/ZP_kn_2.5_large_box/'; % folder = '/misc/gene_results/CBC/128x64x16x24x12/'; % folder = '/misc/gene_results/CBC/196x96x20x32x16_02/'; % folder = '/misc/gene_results/CBC/128x64x16x6x4/'; % folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_01/'; % folder = '/misc/gene_results/CBC/KT_4.5_128x64x16x24x12_01/'; % folder = '/misc/gene_results/CBC/KT_9_128x64x16x24x12/'; % folder = '/misc/gene_results/CBC/KT_13_large_box_128x64x16x24x12/'; % folder = '/misc/gene_results/CBC/Lapillone_Fig6/'; -folder = '/misc/gene_results/Z-pinch/HP_kN_1.6_adapt_mu_01/'; +% folder = '/misc/gene_results/Z-pinch/HP_kN_1.6_adapt_mu_01/'; +folder = '/misc/gene_results/miller/'; gene_data = load_gene_data(folder); gene_data = invert_kxky_to_kykx_gene_results(gene_data); if 1 %% Space time diagramm (fig 11 Ivanov 2020) options.TAVG_0 = 0.1*gene_data.Ts3D(end); options.TAVG_1 = gene_data.Ts3D(end); % Averaging times duration options.NMVA = 1; % Moving average for time traces options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x, Q_x) options.INTERP = 1; gene_data.FIGDIR = folder; fig = plot_radial_transport_and_spacetime(gene_data,options); save_figure(gene_data,fig,'.png') end if 0 %% statistical transport averaging options.T = [100 500]; fig = statistical_transport_averaging(gene_data,options); end if 0 %% 2D snapshots % Options options.INTERP = 1; options.POLARPLOT = 0; options.AXISEQUAL = 0; % options.NAME = 'Q_x'; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME = '\Gamma_x'; % options.NAME = 'k^2n_e'; options.PLAN = 'kxky'; % options.NAME ='f_e'; % options.PLAN = 'sx'; options.COMP = 'avg'; options.TIME = [325 400]; gene_data.a = data.EPS * 2000; fig = photomaton(gene_data,options); save_figure(gene_data,fig,'.png') end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options options.INTERP = 1; options.POLARPLOT = 0; options.NAME = '\phi'; % options.NAME = 'v_y'; % options.NAME = 'n_i^{NZ}'; % options.NAME = '\Gamma_x'; % options.NAME = 'n_i'; options.PLAN = 'xy'; % options.NAME = 'f_e'; % options.PLAN = 'sx'; options.COMP = 'avg'; options.TIME = gene_data.Ts3D; gene_data.a = data.EPS * 2000; create_film(gene_data,options,'.gif') end if 0 %% Geometry names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',... '$B_0$','$\partial_x B_0$','$\partial_y B_0$','$\partial_z B_0$',... '$J$','$R$','$\phi$','$Z$','$\partial_R x$','$\partial_Z x$'}; figure; subplot(311) for i = 1:6 plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; end xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); title('GENE geometry'); subplot(312) for i = 7:10 plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; end xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); subplot(313) for i = 11:16 plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on; end xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); end if 0 %% Show f_i(vpar,mu) options.times = 250:500; options.specie = 'i'; options.PLT_FCT = 'contour'; options.folder = folder; options.iz = 'avg'; options.FIELD = ''; options.ONED = 0; % options.FIELD = 'Q_es'; plot_fa_gene(options); end if 0 %% Time averaged spectrum options.TIME = 300:600; options.NORM =1; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME ='\Gamma_x'; options.PLAN = 'kxky'; options.COMPZ = 'avg'; options.OK = 0; options.COMPXY = 'zero'; options.COMPT = 'avg'; options.PLOT = 'semilogy'; fig = spectrum_1D(gene_data,options); % save_figure(data,fig) end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = 1; options.TIME = 100:700; options.NMA = 1; options.NMODES = 15; options.iz = 'avg'; fig = mode_growth_meter(gene_data,options); save_figure(gene_data,fig) end diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m index eecf714..65ba29b 100644 --- a/wk/analysis_gyacomo.m +++ b/wk/analysis_gyacomo.m @@ -1,222 +1,228 @@ %% UNCOMMENT FOR TUTORIAL % gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory % resdir = '.'; %Name of the directory where the results are located % JOBNUMMIN = 00; JOBNUMMAX = 10; %% addpath(genpath([gyacomodir,'matlab'])) % ... add addpath(genpath([gyacomodir,'matlab/plot'])) % ... add addpath(genpath([gyacomodir,'matlab/compute'])) % ... add addpath(genpath([gyacomodir,'matlab/load'])) % ... add %% Load the results LOCALDIR = [gyacomodir,resdir,'/']; MISCDIR = ['/misc/gyacomo_outputs/',resdir,'/']; %For long term storage system(['mkdir -p ',MISCDIR]); system(['mkdir -p ',LOCALDIR]); CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD); system(CMD); % Load outputs from jobnummin up to jobnummax data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing data.localdir = LOCALDIR; data.FIGDIR = LOCALDIR; %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 1 %% Space time diagramm (fig 11 Ivanov 2020) % data.scale = 1;%/(data.Nx*data.Ny)^2; -i_ = 11; +i_ = 1; disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))]) disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))]) options.TAVG_0 = data.TJOB_SE(i_);%0.4*data.Ts3D(end); options.TAVG_1 = data.TJOB_SE(i_+1);%0.9*data.Ts3D(end); % Averaging times duration options.NCUT = 4; % Number of cuts for averaging and error estimation -options.NMVA = 100; % Moving average for time traces +options.NMVA = 1; % Moving average for time traces % options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) -% options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) -% options.INTERP = 1; +options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) +options.INTERP = 0; +options.RESOLUTION = 256; fig = plot_radial_transport_and_spacetime(data,options); save_figure(data,fig,'.png') end if 0 %% statistical transport averaging -options.T = [200 400]; +for i_ = 1:2:21 +% i_ = 3; +disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))]) +disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))]) +options.T = [data.TJOB_SE(i_) data.TJOB_SE(i_+1)]; +options.NPLOTS = 0; fig = statistical_transport_averaging(data,options); end - +end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options options.INTERP = 1; options.POLARPLOT = 0; options.NAME = '\phi'; % options.NAME = '\omega_z'; -% options.NAME = 'N_i^{00}'; +% options.NAME = 'N_e^{00}'; % options.NAME = 'v_y'; % options.NAME = 'n_i^{NZ}'; % options.NAME = '\Gamma_x'; % options.NAME = 'n_i'; options.PLAN = 'xy'; % options.NAME = 'f_i'; % options.PLAN = 'sx'; -% options.COMP = 'avg'; +options.COMP = 'avg'; % options.TIME = data.Ts5D(end-30:end); % options.TIME = data.Ts3D; -options.TIME = [000:50:7000]; +options.TIME = [000:0.1:7000]; data.EPS = 0.1; data.a = data.EPS * 2000; options.RESOLUTION = 256; create_film(data,options,'.gif') end if 1 %% 2D snapshots % Options options.INTERP = 1; options.POLARPLOT = 0; options.AXISEQUAL = 1; options.NAME = '\phi'; % options.NAME = '\psi'; % options.NAME = 'n_i'; % options.NAME = 'N_i^{00}'; % options.NAME = 'T_i'; % options.NAME = '\Gamma_x'; % options.NAME = 'k^2n_e'; options.PLAN = 'xy'; % options.NAME 'f_i'; % options.PLAN = 'sx'; options.COMP = 'avg'; options.TIME = [1000 1800 2500 3000 4000]; data.a = data.EPS * 2e3; fig = photomaton(data,options); % save_figure(data,fig) end if 0 %% 3D plot on the geometry options.INTERP = 0; options.NAME = '\phi'; options.PLANES = [1]; options.TIME = [30]; options.PLT_MTOPO = 1; options.PLT_FTUBE = 0; data.EPS = 0.4; data.rho_o_R = 3e-3; % Sound larmor radius over Machine size ratio fig = show_geometry(data,options); save_figure(data,fig,'.png') end if 0 %% Kinetic distribution function sqrt(xy) (GENE vsp) options.SPAR = linspace(-3,3,32)+(6/127/2); options.XPERP = linspace( 0,6,32); % options.SPAR = gene_data.vp'; % options.XPERP = gene_data.mu'; options.iz = 'avg'; options.T = [250 600]; options.PLT_FCT = 'pcolor'; options.ONED = 0; options.non_adiab = 0; options.SPECIE = 'i'; options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene fig = plot_fa(data,options); % save_figure(data,fig,'.png') end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 0; options.ST = 1; options.PLOT_TYPE = 'space-time'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [1000]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); % save_figure(data,fig,'.png'); end if 0 %% Time averaged spectrum -options.TIME = [300 600]; +options.TIME = [2000 3000]; options.NORM =1; options.NAME = '\phi'; % options.NAME = 'N_i^{00}'; % options.NAME ='\Gamma_x'; options.PLAN = 'kxky'; options.COMPZ = 'avg'; options.OK = 0; options.COMPXY = 'avg'; % avg/sum/max/zero/ 2D plot otherwise options.COMPT = 'avg'; options.PLOT = 'semilogy'; fig = spectrum_1D(data,options); % save_figure(data,fig,'.png') end if 0 %% 1D real plot options.TIME = [50 100 200]; options.NORM = 0; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME ='\Gamma_x'; % options.NAME ='s_y'; options.COMPX = 'avg'; options.COMPY = 'avg'; options.COMPZ = 1; options.COMPT = 1; options.MOVMT = 1; fig = real_plot_1D(data,options); % save_figure(data,fig,'.png') end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = [0.1 0.2 0.3 0.4]; options.TIME = [00:1200]; options.NMA = 1; options.NMODES = 5; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end if 0 %% ZF caracteristics (space time diagrams) TAVG_0 = 1200; TAVG_1 = 1500; % Averaging times duration % chose your field to plot in spacetime diag (uzf,szf,Gx) fig = ZF_spacetime(data,TAVG_0,TAVG_1); save_figure(data,fig,'.png') end if 0 %% Metric infos fig = plot_metric(data); end if 0 %% linear growth rate for 3D fluxtube trange = [0 100]; nplots = 1; lg = compute_fluxtube_growth_rate(data,trange,nplots); end if 0 %% linear growth rate for 3D Zpinch trange = [5 15]; options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 1; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig,'.png') end diff --git a/wk/header_2DZP_results.m b/wk/header_2DZP_results.m index b180dd5..10428ea 100644 --- a/wk/header_2DZP_results.m +++ b/wk/header_2DZP_results.m @@ -1,200 +1,203 @@ %% Directory of the simulation gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory % if 1% Local results resdir =''; resdir =''; resdir =''; resdir =''; % resdir ='debug/ppj_init'; %% nu = 5e-1 % Sugama % resdir ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_5e-01_SGGK';% also in 7x4 % resdir ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_5e-01_SGGK'; % resdir ='Hallenbert_nu_5e-01/200x32_7x4_L_120_kN_1.7_kT_0.425_nu_5e-01_SGGK';% also in 7x4 % resdir ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_5e-01_SGGK'; % resdir ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_5e-01_SGGK';%also in 7x4 %% nu = 1e-1 % Landau % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_LDGK'; % resdir ='Hallenbert_nu_1e-01/150x50_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_LDGK'; % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_LDGK'; % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_LDGK'; % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_LDGK'; % resdir ='Hallenbert_nu_1e-01/150x50_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_LDGK'; % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-01_LDGK'; % Sugama % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_SGGK'; % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_SGGK'; % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_SGGK'; % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-01_SGGK'; % Dougherty % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_DGGK'; % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_DGGK'; % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_DGGK'; % resdir ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_DGGK'; %% nu = 5e-2 % resdir ='Hallenbert_nu_5e-02/200x32_11x6_L_120_kN_1.8_kT_0.45_nu_5e-02_SGGK';%For GENE benchmark % to analyse (added HD) % resdir ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGGK'; % testing various NL closures % resdir ='Hallenbert_nu_5e-02/200x32_7x4_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGGK'; % resdir ='Hallenbert_nu_5e-02/200x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK'; % resdir ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK'; % resdir ='Hallenbert_nu_5e-02/256x64_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK'; % resdir ='Hallenbert_nu_5e-02/256x64_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK'; % resdir ='Hallenbert_nu_5e-02/200x32_21x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK'; % resdir ='Hallenbert_nu_5e-02/200x32_17x9_L_120_kN_1.8_kT_0.45_nu_5e-02_SGDK'; % resdir ='Hallenbert_nu_5e-02/200x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_DGGK'; % resdir ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_DGGK'; % resdir ='Hallenbert_nu_5e-02/128x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_FCGK'; %% nu = 1e-2 % Landau % resdir ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-02_LDGK'; % resdir ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_LDGK'; % resdir ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_LDGK'; % resdir ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_LDGK'; % resdir ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-02_LDGK'; % Sugama % resdir ='kobayashi_2015_fig1/150x150_5x3_L_100_kN_1.4_nu_5e-03_SGGK'; % resdir ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.5_kT_0.375_nu_1e-02_SGGK'; % resdir ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-02_SGGK'; % resdir ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_SGGK'; % resdir ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.6_kT_0.4_nu_1e-02_SGGK'; % resdir ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK'; % resdir ='Hallenbert_nu_1e-02/300x64_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK'; % resdir ='Hallenbert_nu_1e-02/200x32_11x6_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK'; % resdir ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_SGGK'; % To analyse (added HD) % resdir ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-02_SGGK'; % resdir ='Hallenbert_nu_1e-02/200x32_11x6_L_120_kN_1.9_kT_0.475_nu_1e-02_SGGK'; % Dougherty % resdir ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.5_kT_0.375_nu_1e-02_DGGK'; % resdir ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_DGGK'; % resdir ='Hallenbert_nu_1e-02/200x32_8x5_L_120_kN_1.6_kT_0.4_nu_0e+00_DGGK'; % resdir ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_DGGK'; %% nu = 5e-3 % resdir ='Hallenbert_nu_5e-03/200x32_5x3_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_muxy_5e-2'; % resdir ='Hallenbert_nu_5e-03/200x32_5x3_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_mux_5e-2_muy_6e-1'; % resdir ='Hallenbert_nu_5e-03/200x32_11x6_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_mux_5e-2_muy_6e-1'; % resdir ='Hallenbert_nu_5e-03/200x32_11x6_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_muxy_5e-2'; %% nu = 0 % resdir ='Hallenbert_fig2a/200x32_21x11_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2'; % resdir ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2'; % resdir ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2'; % resdir ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuDGGK_0.1_muxy_1e-2'; % resdir ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nuDGGK_0.1_muxy_1e-2'; % resdir ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuSGGK_0.1_muxy_1e-2'; % resdir ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nuSGGK_0.1_muxy_1e-2'; % resdir ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuLDGK_0.1_muxy_1e-2'; % resdir ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuLRGK_0.1_muxy_1e-2'; % resdir ='Hallenbert_fig2b/200x32_11x6_Lx_240_Ly_120_kN_2.5_eta_0.25_nu_0_muxy_1e-1'; % resdir ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nu_0_muxy_1e-1'; % resdir ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuSGGK_0.1_muxy_1e-1'; % resdir ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuDGGK_0.1_muxy_1e-1'; % resdir ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuLDGK_0.1_muxy_1e-1'; % resdir ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuLRGK_0.1_muxy_1e-1'; % resdir ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nu_0_muxy_5e-2'; % resdir ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nu_0_muxy_5e-2'; % resdir ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuSGGK_0.1_muxy_5e-2'; % resdir ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nuSGGK_0.1_muxy_5e-2'; % resdir ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuDGGK_0.1_muxy_5e-2'; % resdir ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nuDGGK_0.1_muxy_5e-2'; % resdir ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuLDGK_0.1_muxy_5e-2'; % resdir ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuLRGK_0.1_muxy_5e-2'; % resdir ='Hallenbert_fig2c/200x32_9x5_Lx_120_Ly_60_kN_2.0_eta_0.25_nuLRGK_0.1_muxy_5e-2'; %% Transport scan % resdir = 'nu_0.1_transport_scan/colless_kn_1.7_to_2.0'; % resdir = 'nu_0.1_transport_scan/colless_kn_2.1_to_2.5'; % resdir = 'nu_0.1_transport_scan/LB_kn_2.0'; % resdir = 'nu_0.1_transport_scan/DG_kn_1.8_to_2.1'; % resdir = 'nu_0.1_transport_scan/DG_kn_2.2_to_2.5'; % resdir = 'nu_0.1_transport_scan/DG_conv_kN_1.9'; % resdir = 'nu_0.1_transport_scan/SG_kn_1.7_to_2.0'; % resdir = 'nu_0.1_transport_scan/SG_10x5_conv_test'; % resdir = 'nu_0.1_transport_scan/SG_kn_2.2_to_2.5'; % resdir = 'nu_0.1_transport_scan/LD_kn_2.0_to_2.5'; % resdir = 'nu_0.1_transport_scan/LD_kn_1.7_to_2.5'; % resdir = 'nu_0.1_transport_scan/LR_kn_1.7_to_2.0'; % resdir = 'nu_0.1_transport_scan/LR_kn_2.1_to_2.5'; % resdir = 'nu_0.1_transport_scan/colless_kn_2.2_Lx1.5'; % resdir = 'nu_0.1_transport_scan/colless_kn_2.2_HD'; % resdir = 'nu_0.1_transport_scan/colless_kn_1.6_HD'; % resdir = 'nu_0.1_transport_scan/large_box_kN_2.1_nu_0.1'; % resdir = 'nu_0.1_transport_scan/large_box_kN_2.0_nu_0.1'; % resdir = 'predator_prey_nu_scan/DG_Kn_1.7_nu_0.01'; % resdir = 'ZF_damping_linear_nu_0_20x10_kn_1.6_GK/LR_4x2_nu_0.1'; % resdir = 'ZF_damping_nu_0_20x10_kn_1.6_GK/HSG_4x2_nu_0.1'; % resdir = 'ZF_damping_nu_0_5x3_kn_2.5_GK/LR_4x2_nu_0.1'; % resdir = 'hacked_sugama/hacked_B_kn_1.6_200x32_L_120x60_nu_0.1'; % resdir = 'shearless_cyclone/200x32x24_5x4_Lx_120_Ly_60_q0_1.4_e_0.18_kN_2.22_kT_6.9_nuLR_0.01_adiab_e'; % resdir = 'shearless_cyclone/no_sg_128x32x36_6x3_Lx_120_Ly_60_q0_1.4_e_0.18_kN_2.22_kT_6.9_adiab_e'; % resdir = 'shearless_cyclone/sgrid_128x64x32_4x2_Lx_100_Ly_120_q0_1.4_e_0.18_kN_2.22_kT_6.96_adiab_e'; % resdir = 'shearless_cyclone/sgrid_128x64x32_4x2_Lx_100_Ly_120_q0_1.4_e_0.18_kN_1.78_kT_5.52_adiab_e'; % resdir = 'linear_shearless_cyclone/4_2_cyclone_1.0'; % resdir = 'linear_shearless_cyclone/test_fmom'; % else% Marconi results % resdir =''; % resdir =''; % resdir ='';fd % resdir =''; % resdir ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A_new/300x300_5x3_L_120_kN_1.6667_nu_1e-01_SGGK/out.txt'; % % resdir ='/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'; % % BASIC.RESDIR = ['../',resdir(46:end-8),'/']; % MISCDIR = ['/misc/HeLaZ_outputs/',resdir(46:end-8),'/']; % end %% ZPINCH rerun % resdir ='Zpinch_rerun/Kn_2.5_200x48x5x3'; % resdir ='Zpinch_rerun/Kn_2.5_256x128x5x3'; % resdir ='Zpinch_rerun/Kn_2.5_312x196x5x3_Lx_400_Ly_200'; % resdir ='Zpinch_rerun/Kn_2.5_256x64x5x3'; % resdir ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box'; % resdir ='Zpinch_rerun/Kn_2.0_256x64x9x5_Lx_240_Ly_120'; % resdir ='Zpinch_rerun/Kn_1.6_256x128x7x4'; % resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6'; % resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6_conv'; % resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6_mu_0.5'; % resdir ='Zpinch_rerun/Kn_1.6_256x128x21x11'; %% nu scan % resdir = 'Zpinch_rerun/kN_2.2_coll_scan_128x48x5x3'; % resdir = 'Zpinch_rerun/Ultra_HD_312x196x5x3'; % resdir = 'Zpinch_rerun/UHD_nu_001_LDGK'; % resdir = 'Zpinch_rerun/UHD_nu_01_LDGK'; % resdir = 'Zpinch_rerun/UHD_nu_1_LDGK'; % resdir ='Zpinch_rerun/kN_1.7_SGGK_conv_200x32x7x3_nu_0.01'; % resdir ='Zpinch_rerun/kN_1.7_LDGKii_200x32x7x3_nu_scan'; % resdir ='Zpinch_rerun/nu_0.1_LDGKii_200x48x7x4_kN_scan'; -resdir ='Zpinch_rerun/nu_0.1_FCGK_200x48x5x3_kN_scan'; +% resdir ='Zpinch_rerun/nu_0.1_FCGK_200x48x5x3_kN_scan'; +% resdir ='Zpinch_rerun/nu_0.1_SGGK_200x48x5x3_kN_scan'; % resdir = 'Zpinch_rerun/kN_1.7_FCGK_200x32x5x3_nu_scan'; % resdir = 'Zpinch_rerun/kN_1.7_SGGK_200x32x7x4_nu_scan'; +% resdir = 'Zpinch_rerun/kN_1.7_SGGK_200x32x5x3_nu_scan'; +resdir = 'Zpinch_rerun/kN_1.7_SGGK_256x64x5x3_nu_scan'; %% JOBNUMMIN = 00; JOBNUMMAX = 10; resdir = ['results/',resdir]; run analysis_gyacomo \ No newline at end of file diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m index 2a56c2d..189cc0e 100644 --- a/wk/header_3D_results.m +++ b/wk/header_3D_results.m @@ -1,64 +1,64 @@ % Directory of the code "mypathtoHeLaZ/HeLaZ/" -helazdir = '/home/ahoffman/HeLaZ/'; +gyacomodir = '/home/ahoffman/gyacomo/'; % Directory of the simulation (from results) % if 1% Local results -% outfile ='volcokas/64x32x16x5x3_kin_e_npol_1'; +% resdir ='volcokas/64x32x16x5x3_kin_e_npol_1'; %% Dimits -% outfile ='shearless_cyclone/128x64x16x5x3_Dim_90'; -% outfile ='shearless_cyclone/128x64x16x9x5_Dim_scan/128x64x16x9x5_Dim_60'; -% outfile ='shearless_cyclone/128x64x16x5x3_Dim_scan/128x64x16x5x3_Dim_70'; -% outfile ='shearless_cyclone/64x32x16x5x3_Dim_scan/64x32x16x5x3_Dim_70'; +% resdir ='shearless_cyclone/128x64x16x5x3_Dim_90'; +% resdir ='shearless_cyclone/128x64x16x9x5_Dim_scan/128x64x16x9x5_Dim_60'; +% resdir ='shearless_cyclone/128x64x16x5x3_Dim_scan/128x64x16x5x3_Dim_70'; +% resdir ='shearless_cyclone/64x32x16x5x3_Dim_scan/64x32x16x5x3_Dim_70'; %% AVS -% outfile = 'volcokas/64x32x16x5x3_kin_e_npol_1'; -% outfile = 'volcokas/64x32x16x5x3_kin_e_npol_1'; -% outfile = 'shearless_cyclone/64x32x80x5x3_CBC_Npol_5_kine'; -% outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; -% outfile = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine'; -% outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; +% resdir = 'volcokas/64x32x16x5x3_kin_e_npol_1'; +% resdir = 'volcokas/64x32x16x5x3_kin_e_npol_1'; +% resdir = 'shearless_cyclone/64x32x80x5x3_CBC_Npol_5_kine'; +% resdir = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; +% resdir = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine'; +% resdir = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; %% shearless CBC -% outfile ='shearless_cyclone/64x32x16x5x3_CBC_080'; -% outfile ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100'; -% outfile ='shearless_cyclone/64x32x16x5x3_CBC_120'; +% resdir ='shearless_cyclone/64x32x16x5x3_CBC_080'; +% resdir ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100'; +% resdir ='shearless_cyclone/64x32x16x5x3_CBC_120'; -% outfile ='shearless_cyclone/64x32x16x9x5_CBC_080'; -% outfile ='shearless_cyclone/64x32x16x9x5_CBC_100'; -% outfile ='shearless_cyclone/64x32x16x9x5_CBC_120'; +% resdir ='shearless_cyclone/64x32x16x9x5_CBC_080'; +% resdir ='shearless_cyclone/64x32x16x9x5_CBC_100'; +% resdir ='shearless_cyclone/64x32x16x9x5_CBC_120'; -% outfile = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK'; +% resdir = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK'; %% CBC -% outfile = 'CBC/64x32x16x5x3'; -% outfile = 'CBC/64x128x16x5x3'; -% outfile = 'CBC/128x64x16x5x3'; -% outfile = 'CBC/96x96x16x3x2_Nexc_6'; -% outfile = 'CBC/128x96x16x3x2'; -% outfile = 'CBC/192x96x16x3x2'; -% outfile = 'CBC/192x96x24x13x7'; -% outfile = 'CBC/kT_11_128x64x16x5x3'; -% outfile = 'CBC/kT_9_256x128x16x3x2'; -% outfile = 'CBC/kT_4.5_128x64x16x13x3'; -% outfile = 'CBC/kT_4.5_192x96x24x13x7'; -% outfile = 'CBC/kT_4.5_128x64x16x13x7'; -% outfile = 'CBC/kT_4.5_128x96x24x15x5'; -% outfile = 'CBC/kT_5.3_192x96x24x13x7'; -% outfile = 'CBC/kT_13_large_box_128x64x16x5x3'; -% outfile = 'CBC/kT_11_96x64x16x5x3_ky_0.02'; +% resdir = 'CBC/64x32x16x5x3'; +% resdir = 'CBC/64x128x16x5x3'; +% resdir = 'CBC/128x64x16x5x3'; +% resdir = 'CBC/96x96x16x3x2_Nexc_6'; +% resdir = 'CBC/128x96x16x3x2'; +% resdir = 'CBC/192x96x16x3x2'; +% resdir = 'CBC/192x96x24x13x7'; +% resdir = 'CBC/kT_11_128x64x16x5x3'; +% resdir = 'CBC/kT_9_256x128x16x3x2'; +% resdir = 'CBC/kT_4.5_128x64x16x13x3'; +% resdir = 'CBC/kT_4.5_192x96x24x13x7'; +% resdir = 'CBC/kT_4.5_128x64x16x13x7'; +% resdir = 'CBC/kT_4.5_128x96x24x15x5'; +% resdir = 'CBC/kT_5.3_192x96x24x13x7'; +% resdir = 'CBC/kT_13_large_box_128x64x16x5x3'; +% resdir = 'CBC/kT_11_96x64x16x5x3_ky_0.02'; -% outfile = 'CBC/kT_scan_128x64x16x5x3'; -% outfile = 'CBC/kT_scan_192x96x16x3x2'; -% outfile = 'CBC/kT_13_96x96x16x3x2_Nexc_6'; -% outfile = 'dbg/nexc_dbg'; -outfile = 'CBC/NM_F4_kT_4.5_192x64x24x6x4'; +% resdir = 'CBC/kT_scan_128x64x16x5x3'; +% resdir = 'CBC/kT_scan_192x96x16x3x2'; +% resdir = 'CBC/kT_13_96x96x16x3x2_Nexc_6'; +% resdir = 'dbg/nexc_dbg'; +% resdir = 'CBC/NM_F4_kT_4.5_192x64x24x6x4'; -% outfile = 'CBC_Ke_EM/192x96x24x5x3'; -% outfile = 'CBC_Ke_EM/96x48x16x5x3'; -% outfile = 'CBC_Ke_EM/minimal_res'; +% resdir = 'CBC_Ke_EM/192x96x24x5x3'; +% resdir = 'CBC_Ke_EM/96x48x16x5x3'; +% resdir = 'CBC_Ke_EM/minimal_res'; %% KBM -% outfile = 'NL_KBM/192x64x24x5x3'; +% resdir = 'NL_KBM/192x64x24x5x3'; %% Linear CBC -% outfile = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe'; - +% resdir = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe'; +resdir = 'testcases/miller_example_bug'; JOBNUMMIN = 00; JOBNUMMAX = 10; -run analysis_HeLaZ +run analysis_gyacomo diff --git a/wk/lin_ITG.m b/wk/lin_EPY.m similarity index 74% copy from wk/lin_ITG.m copy to wk/lin_EPY.m index bdf3ec0..8a00601 100644 --- a/wk/lin_ITG.m +++ b/wk/lin_EPY.m @@ -1,186 +1,178 @@ -%% QUICK RUN SCRIPT +%% QUICK RUN SCRIPT for linear entropy mode in a Zpinch % This script create a directory in /results and run a simulation directly % from matlab framework. It is meant to run only small problems in linear % for benchmark and debugging purpose since it makes matlab "busy" % -SIMID = 'lin_ITG'; % Name of the simulation +SIMID = 'lin_EPY'; % Name of the simulation RUN = 1; % To run or just to load addpath(genpath('../matlab')) % ... add default_plots_options -HELAZDIR = '/home/ahoffman/HeLaZ/'; -% EXECNAME = 'helaz3'; -EXECNAME = 'helaz3_dbg'; +PROGDIR = '/home/ahoffman/gyacomo/'; +EXECNAME = 'gyacomo'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS -NU = 0.0; % Collision frequency +NU = 0.1; % Collision frequency TAU = 1.0; % e/i temperature ratio -K_Ne = 2.22; % ele Density ''' -K_Te = 6.96; % ele Temperature ''' -K_Ni = 2.22; % ion Density gradient drive -K_Ti = 6.96; % ion Temperature ''' -SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) -% SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) -KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons +K_Ne = 2.0; % ele Density ''' +K_Te = 0.5; % ele Temperature ''' +K_Ni = 2.0; % ion Density gradient drive +K_Ti = 0.5; % ion Temperature ''' +SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) +KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons BETA = 0.0; % electron plasma beta %% GRID PARAMETERS P = 4; J = P/2; PMAXE = P; % Hermite basis size of electrons JMAXE = J; % Laguerre " PMAXI = P; % " ions JMAXI = J; % " -NX = 11; % real space x-gridpoints -NY = 2; % '' y-gridpoints +NX = 2; % real space x-gridpoints +NY = 100; % '' y-gridpoints LX = 2*pi/0.8; % Size of the squared frequency domain -LY = 2*pi/0.3; % Size of the squared frequency domain -NZ = 16; % number of perpendicular planes (parallel grid) +LY = 2*pi/0.05; % Size of the squared frequency domain +NZ = 1; % number of perpendicular planes (parallel grid) NPOL = 1; SG = 0; % Staggered z grids option %% GEOMETRY -% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps -GEOMETRY= 's-alpha'; -% GEOMETRY= 'circular'; +GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps Q0 = 1.4; % safety factor SHEAR = 0.8; % magnetic shear NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) EPS = 0.18; % inverse aspect ratio %% TIME PARMETERS -TMAX = 25; % Maximal time unit -DT = 5e-3; % Time step +TMAX = 50; % Maximal time unit +DT = 1e-2; % Time step SPS0D = 1; % Sampling per time unit for 2D arrays SPS2D = 0; % Sampling per time unit for 2D arrays -SPS3D = 5; % Sampling per time unit for 2D arrays +SPS3D = 2; % Sampling per time unit for 2D arrays SPS5D = 1/5; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints JOB2LOAD= -1; %% OPTIONS LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) % Collision operator % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) -CO = 'DG'; -GKCO = 0; % gyrokinetic operator +CO = 'LD'; +GKCO = 1; % gyrokinetic operator ABCO = 1; % interspecies collisions INIT_ZF = 0; ZF_AMP = 0.0; CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS) KERN = 0; % Kernel model (0 : GK) INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom %% OUTPUTS 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 HD_CO = 0.0; % Hyper diffusivity cutoff ratio MU = 0.0; % Hyperdiffusivity coefficient INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; MU_X = MU; % MU_Y = MU; % N_HD = 4; MU_Z = 2.0; % MU_P = 0.0; % MU_J = 0.0; % LAMBDAD = 0.0; NOISE0 = 0.0e-5; % Init noise amplitude BCKGD0 = 1.0; % Init background GRADB = 1.0; CURVB = 1.0; %%------------------------------------------------------------------------- %% RUN setup % system(['rm fort*.90']); % Run linear simulation if RUN -% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) -% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) -% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk']) - system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk']) -% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk']) + system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk']) end %% Load results %% filename = [SIMID,'/',PARAMS,'/']; -LOCALDIR = [HELAZDIR,'results/',filename,'/']; +LOCALDIR = [PROGDIR,'results/',filename,'/']; % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 00; data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing %% Short analysis if 1 %% linear growth rate (adapted for 2D zpinch and fluxtube) trange = [0.5 1]*data.Ts3D(end); -nplots = 3; +nplots = 1; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z lg = compute_fluxtube_growth_rate(data,trange,nplots); [gmax, kmax] = max(lg.g_ky(:,end)); [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky); msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg); msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg); end -if 1 +if 0 %% Ballooning plot options.time_2_plot = [120]; options.kymodes = [0.9]; options.normalized = 1; % options.field = 'phi'; fig = plot_ballooning(data,options); end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 1; options.ST = 1; options.PLOT_TYPE = 'space-time'; % options.PLOT_TYPE = 'Tavg-1D'; % options.PLOT_TYPE = 'Tavg-2D'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [0 50]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); save_figure(data,fig) end if 0 %% linear growth rate for 3D Zpinch (kz fourier transform) trange = [0.5 1]*data.Ts3D(end); options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 0; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig) end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = 1; options.TIME = [0:1000]; options.NMA = 1; options.NMODES = 1; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end -if 1 +if 0 %% RH TEST ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end); [~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D)); plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3)); figure plot(data.Ts3D(it0:it1), plt(data.PHI)); xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$') title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky))) end diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m index bdf3ec0..cffb9cc 100644 --- a/wk/lin_ITG.m +++ b/wk/lin_ITG.m @@ -1,186 +1,184 @@ %% QUICK RUN SCRIPT % This script create a directory in /results and run a simulation directly % from matlab framework. It is meant to run only small problems in linear % for benchmark and debugging purpose since it makes matlab "busy" % SIMID = 'lin_ITG'; % Name of the simulation RUN = 1; % To run or just to load addpath(genpath('../matlab')) % ... add default_plots_options -HELAZDIR = '/home/ahoffman/HeLaZ/'; -% EXECNAME = 'helaz3'; -EXECNAME = 'helaz3_dbg'; +HELAZDIR = '/home/ahoffman/gyacomo/'; +EXECNAME = 'gyacomo'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 0.0; % Collision frequency TAU = 1.0; % e/i temperature ratio K_Ne = 2.22; % ele Density ''' K_Te = 6.96; % ele Temperature ''' K_Ni = 2.22; % ion Density gradient drive K_Ti = 6.96; % ion Temperature ''' SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) % SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons BETA = 0.0; % electron plasma beta %% GRID PARAMETERS P = 4; J = P/2; PMAXE = P; % Hermite basis size of electrons JMAXE = J; % Laguerre " PMAXI = P; % " ions JMAXI = J; % " NX = 11; % real space x-gridpoints -NY = 2; % '' y-gridpoints +NY = 32; % '' y-gridpoints LX = 2*pi/0.8; % Size of the squared frequency domain -LY = 2*pi/0.3; % Size of the squared frequency domain +LY = 2*pi/0.05; % Size of the squared frequency domain NZ = 16; % number of perpendicular planes (parallel grid) NPOL = 1; SG = 0; % Staggered z grids option %% GEOMETRY -% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps GEOMETRY= 's-alpha'; -% GEOMETRY= 'circular'; +% GEOMETRY= 'miller'; Q0 = 1.4; % safety factor SHEAR = 0.8; % magnetic shear NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) EPS = 0.18; % inverse aspect ratio %% TIME PARMETERS -TMAX = 25; % Maximal time unit -DT = 5e-3; % Time step +TMAX = 50; % Maximal time unit +DT = 1e-2; % Time step SPS0D = 1; % Sampling per time unit for 2D arrays SPS2D = 0; % Sampling per time unit for 2D arrays SPS3D = 5; % Sampling per time unit for 2D arrays SPS5D = 1/5; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints JOB2LOAD= -1; %% OPTIONS LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) % Collision operator % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) CO = 'DG'; GKCO = 0; % gyrokinetic operator ABCO = 1; % interspecies collisions INIT_ZF = 0; ZF_AMP = 0.0; CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS) KERN = 0; % Kernel model (0 : GK) INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom %% OUTPUTS 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 HD_CO = 0.0; % Hyper diffusivity cutoff ratio MU = 0.0; % Hyperdiffusivity coefficient INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; MU_X = MU; % MU_Y = MU; % N_HD = 4; MU_Z = 2.0; % MU_P = 0.0; % MU_J = 0.0; % LAMBDAD = 0.0; NOISE0 = 0.0e-5; % Init noise amplitude BCKGD0 = 1.0; % Init background GRADB = 1.0; CURVB = 1.0; %%------------------------------------------------------------------------- %% RUN setup % system(['rm fort*.90']); % Run linear simulation if RUN % system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk']) system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk']) end %% Load results %% filename = [SIMID,'/',PARAMS,'/']; LOCALDIR = [HELAZDIR,'results/',filename,'/']; % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 00; data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing %% Short analysis if 1 %% linear growth rate (adapted for 2D zpinch and fluxtube) trange = [0.5 1]*data.Ts3D(end); nplots = 3; lg = compute_fluxtube_growth_rate(data,trange,nplots); [gmax, kmax] = max(lg.g_ky(:,end)); [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky); msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg); msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg); end if 1 %% Ballooning plot options.time_2_plot = [120]; options.kymodes = [0.9]; options.normalized = 1; % options.field = 'phi'; fig = plot_ballooning(data,options); end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 1; options.ST = 1; options.PLOT_TYPE = 'space-time'; % options.PLOT_TYPE = 'Tavg-1D'; % options.PLOT_TYPE = 'Tavg-2D'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [0 50]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); save_figure(data,fig) end if 0 %% linear growth rate for 3D Zpinch (kz fourier transform) trange = [0.5 1]*data.Ts3D(end); options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 0; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig) end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = 1; options.TIME = [0:1000]; options.NMA = 1; options.NMODES = 1; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end if 1 %% RH TEST ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end); [~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D)); plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3)); figure plot(data.Ts3D(it0:it1), plt(data.PHI)); xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$') title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky))) end