diff --git a/testcases/ExB_shear_demo/2D_Z_pinch/GS2_method/fort_01.90 b/demo/ExB_shear/2D_Z_pinch/GS2_method/fort_01.90 similarity index 100% rename from testcases/ExB_shear_demo/2D_Z_pinch/GS2_method/fort_01.90 rename to demo/ExB_shear/2D_Z_pinch/GS2_method/fort_01.90 diff --git a/demo/ExB_shear/2D_Z_pinch/background_Er/fort_01.90 b/demo/ExB_shear/2D_Z_pinch/background_Er/fort_01.90 new file mode 100644 index 0000000000000000000000000000000000000000..6c7afdd4b280af81cb23e0a4b113a172dd5b6031 --- /dev/null +++ b/demo/ExB_shear/2D_Z_pinch/background_Er/fort_01.90 @@ -0,0 +1,92 @@ +&BASIC + nrun = 1e6 + dt = 0.025 + tmax = 200 + maxruntime = 72000 + job2load = 0 +/ +&GRID + pmax = 2 + jmax = 1 + Nx = 128 + Lx = 200 + Ny = 48 + Ly = 60 + Nz = 1 + SG = .f. + Nexc = 1 +/ +&GEOMETRY + geom = 'z-pinch' +/ +&DIAGNOSTICS + dtsave_0d = 1 + dtsave_1d = -1 + dtsave_2d = 0.5 + dtsave_3d = 0.5 + dtsave_5d = 100 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'nonlinear' + Na = 2 ! number of species + mu_x = 1.0 + mu_y = 1.0 + N_HD = 4 + mu_z = 0.0 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.1 + beta = 0.0 + ADIAB_E = .f. + tau_i = 1.0 + ZFrate = 1.0 + ikxZF = 1 +/ +&CLOSURE + hierarchy_closure='truncation' + !hierarchy_closure='max_degree' + dmax = 2 + nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing) + nmax = -1 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 2.0 + k_T_ = 0.4 +/ +&SPECIES + ! electrons + name_ = 'electrons' + tau_ = 1.0 + sigma_= 0.023338 + q_ =-1.0 + k_N_ = 2.0 + k_T_ = 0.4 +/ +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .t. +/ +&INITIAL + INIT_OPT = 'phi' !(phi,blob) + init_background = 0.0 + init_noiselvl = 0.005 + iseed = 42 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' +/ diff --git a/demo/ExB_shear/2D_Z_pinch/control/fort_01.90 b/demo/ExB_shear/2D_Z_pinch/control/fort_01.90 new file mode 100644 index 0000000000000000000000000000000000000000..5899fc021b571093f2888bc11c17d898f07800cf --- /dev/null +++ b/demo/ExB_shear/2D_Z_pinch/control/fort_01.90 @@ -0,0 +1,91 @@ +&BASIC + nrun = 1e6 + dt = 0.05 + tmax = 200 + maxruntime = 72000 + job2load = 0 +/ +&GRID + pmax = 2 + jmax = 1 + Nx = 128 + Lx = 200 + Ny = 48 + Ly = 60 + Nz = 1 + SG = .f. + Nexc = 1 +/ +&GEOMETRY + geom = 'z-pinch' +/ +&DIAGNOSTICS + dtsave_0d = 1 + dtsave_1d = -1 + dtsave_2d = 0.5 + dtsave_3d = 0.5 + dtsave_5d = 100 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'nonlinear' + Na = 2 ! number of species + mu_x = 1.0 + mu_y = 1.0 + N_HD = 4 + mu_z = 0.0 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.1 + beta = 0.0 + ADIAB_E = .f. + tau_i = 1.0 + ExBrate = 0 +/ +&CLOSURE + hierarchy_closure='truncation' + !hierarchy_closure='max_degree' + dmax = 2 + nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing) + nmax = -1 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 2.0 + k_T_ = 0.4 +/ +&SPECIES + ! electrons + name_ = 'electrons' + tau_ = 1.0 + sigma_= 0.023338 + q_ =-1.0 + k_N_ = 2.0 + k_T_ = 0.4 +/ +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .t. +/ +&INITIAL + INIT_OPT = 'phi' !(phi,blob) + init_background = 0.0 + init_noiselvl = 0.005 + iseed = 42 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' +/ diff --git a/testcases/ExB_shear_demo/2D_Z_pinch/fort_00.90 b/demo/ExB_shear/2D_Z_pinch/fort_00.90 similarity index 100% rename from testcases/ExB_shear_demo/2D_Z_pinch/fort_00.90 rename to demo/ExB_shear/2D_Z_pinch/fort_00.90 diff --git a/testcases/ExB_shear_demo/GS2_method/fast_analysis.m b/demo/ExB_shear/GS2_method/fast_analysis.m similarity index 100% rename from testcases/ExB_shear_demo/GS2_method/fast_analysis.m rename to demo/ExB_shear/GS2_method/fast_analysis.m diff --git a/testcases/ExB_shear_demo/GS2_method/fort_00.90 b/demo/ExB_shear/GS2_method/fort_00.90 similarity index 100% rename from testcases/ExB_shear_demo/GS2_method/fort_00.90 rename to demo/ExB_shear/GS2_method/fort_00.90 diff --git a/testcases/ExB_shear_demo/background_Er/fast_analysis.m b/demo/ExB_shear/background_Er/fast_analysis.m similarity index 100% rename from testcases/ExB_shear_demo/background_Er/fast_analysis.m rename to demo/ExB_shear/background_Er/fast_analysis.m diff --git a/testcases/ExB_shear_demo/background_Er/fort.90 b/demo/ExB_shear/background_Er/fort.90 similarity index 100% rename from testcases/ExB_shear_demo/background_Er/fort.90 rename to demo/ExB_shear/background_Er/fort.90 diff --git a/testcases/ExB_shear_demo/background_Er/fort_00.90 b/demo/ExB_shear/background_Er/fort_00.90 similarity index 100% rename from testcases/ExB_shear_demo/background_Er/fort_00.90 rename to demo/ExB_shear/background_Er/fort_00.90 diff --git a/demo/Laguerre_aliasing/2D_Z_pinch/Nmax_0/fort_00.90 b/demo/Laguerre_aliasing/2D_Z_pinch/Nmax_0/fort_00.90 new file mode 100644 index 0000000000000000000000000000000000000000..495787124f3c39313bea9ccf7d5566625236a854 --- /dev/null +++ b/demo/Laguerre_aliasing/2D_Z_pinch/Nmax_0/fort_00.90 @@ -0,0 +1,91 @@ +&BASIC + nrun = 1e6 + dt = 0.02 + tmax = 500 + maxruntime = 72000 + job2load = -1 +/ +&GRID + pmax = 2 + jmax = 1 + Nx = 128 + Lx = 200 + Ny = 48 + Ly = 60 + Nz = 1 + SG = .f. + Nexc = 1 +/ +&GEOMETRY + geom = 'z-pinch' +/ +&DIAGNOSTICS + dtsave_0d = 1 + dtsave_1d = -1 + dtsave_2d = -1 + dtsave_3d = 2 + dtsave_5d = 50 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'nonlinear' + Na = 2 ! number of species + mu_x = 1.0 + mu_y = 0.0 + N_HD = 4 + mu_z = 0.0 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.05 + beta = 0.0 + ADIAB_E = .f. + tau_i = 1.0 +/ +&CLOSURE + !hierarchy_closure='truncation' + hierarchy_closure='max_degree' + dmax = 2 + !nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing) + nonlinear_closure='truncation' + nmax = 0 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 2.0 + k_T_ = 0.4 +/ +&SPECIES + ! electrons + name_ = 'electrons' + tau_ = 1.0 + sigma_= 0.023338 + q_ =-1.0 + k_N_ = 2.0 + k_T_ = 0.4 +/ +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .t. +/ +&INITIAL + INIT_OPT = 'phi' !(phi,blob) + init_background = 0.0 + init_noiselvl = 0.005 + iseed = 42 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' +/ diff --git a/demo/Laguerre_aliasing/2D_Z_pinch/anti_Laguerre_aliasing/fort_00.90 b/demo/Laguerre_aliasing/2D_Z_pinch/anti_Laguerre_aliasing/fort_00.90 new file mode 100644 index 0000000000000000000000000000000000000000..79b5231e35d9a595a683c1c4c05f135b43e7edd7 --- /dev/null +++ b/demo/Laguerre_aliasing/2D_Z_pinch/anti_Laguerre_aliasing/fort_00.90 @@ -0,0 +1,90 @@ +&BASIC + nrun = 1e6 + dt = 0.02 + tmax = 500 + maxruntime = 72000 + job2load = -1 +/ +&GRID + pmax = 2 + jmax = 1 + Nx = 128 + Lx = 200 + Ny = 48 + Ly = 60 + Nz = 1 + SG = .f. + Nexc = 1 +/ +&GEOMETRY + geom = 'z-pinch' +/ +&DIAGNOSTICS + dtsave_0d = 1 + dtsave_1d = -1 + dtsave_2d = -1 + dtsave_3d = 2 + dtsave_5d = 50 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'nonlinear' + Na = 2 ! number of species + mu_x = 1.0 + mu_y = 0.0 + N_HD = 4 + mu_z = 0.0 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.05 + beta = 0.0 + ADIAB_E = .f. + tau_i = 1.0 +/ +&CLOSURE + !hierarchy_closure='truncation' + hierarchy_closure='max_degree' + dmax = 2 + nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing) + nmax = -1 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 2.0 + k_T_ = 0.4 +/ +&SPECIES + ! electrons + name_ = 'electrons' + tau_ = 1.0 + sigma_= 0.023338 + q_ =-1.0 + k_N_ = 2.0 + k_T_ = 0.4 +/ +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .t. +/ +&INITIAL + INIT_OPT = 'phi' !(phi,blob) + init_background = 0.0 + init_noiselvl = 0.005 + iseed = 42 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' +/ diff --git a/demo/Laguerre_aliasing/2D_Z_pinch/fort_00.90 b/demo/Laguerre_aliasing/2D_Z_pinch/fort_00.90 new file mode 100644 index 0000000000000000000000000000000000000000..190c9ffbe3b05f2fbce4bc0bab4d2cb822397499 --- /dev/null +++ b/demo/Laguerre_aliasing/2D_Z_pinch/fort_00.90 @@ -0,0 +1,90 @@ +&BASIC + nrun = 1e6 + dt = 0.05 + tmax = 500 + maxruntime = 72000 + job2load = -1 +/ +&GRID + pmax = 2 + jmax = 1 + Nx = 128 + Lx = 200 + Ny = 48 + Ly = 60 + Nz = 1 + SG = .f. + Nexc = 1 +/ +&GEOMETRY + geom = 'z-pinch' +/ +&DIAGNOSTICS + dtsave_0d = 1 + dtsave_1d = -1 + dtsave_2d = 0.5 + dtsave_3d = 0.5 + dtsave_5d = 100 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'nonlinear' + Na = 2 ! number of species + mu_x = 1.0 + mu_y = 0.0 + N_HD = 4 + mu_z = 0.0 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.05 + beta = 0.0 + ADIAB_E = .f. + tau_i = 1.0 +/ +&CLOSURE + hierarchy_closure='truncation' + !hierarchy_closure='max_degree' + dmax = 2 + nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing) + nmax = -1 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 2.0 + k_T_ = 0.4 +/ +&SPECIES + ! electrons + name_ = 'electrons' + tau_ = 1.0 + sigma_= 0.023338 + q_ =-1.0 + k_N_ = 2.0 + k_T_ = 0.4 +/ +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .t. +/ +&INITIAL + INIT_OPT = 'phi' !(phi,blob) + init_background = 0.0 + init_noiselvl = 0.005 + iseed = 42 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' +/ diff --git a/demo/Laguerre_aliasing/2D_Z_pinch/maximized_sum/fort_00.90 b/demo/Laguerre_aliasing/2D_Z_pinch/maximized_sum/fort_00.90 new file mode 100644 index 0000000000000000000000000000000000000000..a16308273acbe715abde82e3ed416bcc65787632 --- /dev/null +++ b/demo/Laguerre_aliasing/2D_Z_pinch/maximized_sum/fort_00.90 @@ -0,0 +1,91 @@ +&BASIC + nrun = 1e6 + dt = 0.02 + tmax = 500 + maxruntime = 72000 + job2load = -1 +/ +&GRID + pmax = 2 + jmax = 1 + Nx = 128 + Lx = 200 + Ny = 48 + Ly = 60 + Nz = 1 + SG = .f. + Nexc = 1 +/ +&GEOMETRY + geom = 'z-pinch' +/ +&DIAGNOSTICS + dtsave_0d = 1 + dtsave_1d = -1 + dtsave_2d = -1 + dtsave_3d = 2 + dtsave_5d = 50 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'nonlinear' + Na = 2 ! number of species + mu_x = 1.0 + mu_y = 0.0 + N_HD = 4 + mu_z = 0.0 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.05 + beta = 0.0 + ADIAB_E = .f. + tau_i = 1.0 +/ +&CLOSURE + !hierarchy_closure='truncation' + hierarchy_closure='max_degree' + dmax = 2 + !nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing) + nonlinear_closure='full_sum' + nmax = -1 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 2.0 + k_T_ = 0.4 +/ +&SPECIES + ! electrons + name_ = 'electrons' + tau_ = 1.0 + sigma_= 0.023338 + q_ =-1.0 + k_N_ = 2.0 + k_T_ = 0.4 +/ +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .t. +/ +&INITIAL + INIT_OPT = 'phi' !(phi,blob) + init_background = 0.0 + init_noiselvl = 0.005 + iseed = 42 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' +/ diff --git a/matlab/compute/compute_fa_2D.m b/matlab/compute/compute_fa_2D.m index 8d9758e0b27a106faad7bce9a9bce4058ee67f57..f2e9552e97653b21476b4888f0fd77b7aaa31f4f 100644 --- a/matlab/compute/compute_fa_2D.m +++ b/matlab/compute/compute_fa_2D.m @@ -1,4 +1,4 @@ -function [SS,XX,FF] = compute_fa_2D(data, options) +function [SS,XX,FF] = compute_fa_2D(data, species, s, x, T) %% Compute the dispersion function from the moment hierarchi decomp. % Normalized Hermite Hp = @(p,s) polyval(HermitePoly(p),s)./sqrt(2.^p.*factorial(p)); @@ -9,99 +9,90 @@ Lj = @(j,x) polyval(LaguerrePoly(j),x); FaM = @(s,x) exp(-s.^2-x); %meshgrid for efficient evaluation -[SS, XX] = meshgrid(options.SPAR, options.XPERP); +[SS, XX] = meshgrid(s, x); -[~, ikx0] = min(abs(data.kx)); -[~, iky0] = min(abs(data.ky)); -kx_ = data.kx; ky_ = data.ky; - -switch options.SPECIES +switch species case 'e' - Napj_ = data.Nepj; - parray = double(data.Pe); - jarray = double(data.Je); - case 'i' - Napj_ = data.Nipj; - parray = double(data.Pi); - jarray = double(data.Ji); -end + Napj_ = data.Napjz(2,:,:,:,:); -switch options.iz - case 'avg' - Napj_ = mean(Napj_,5); - phi_ = mean(data.PHI,3); - otherwise - iz = options.iz; - Napj_ = Napj_(:,:,:,:,iz,:); - phi_ = data.PHI(:,:,iz); + case 'i' + Napj_ = data.Napjz(1,:,:,:,:); end +parray = double(data.grids.Parray); +jarray = double(data.grids.Jarray); +% switch options.iz + % case 'avg' + options.SHOW_FLUXSURF = 0; + options.SHOW_METRICS = 0; + options.SHOW_CURVOP = 0; + [~, geo_arrays] = plot_metric(data,options); + J = geo_arrays.Jacobian; + Nz = data.grids.Nz; + tmp_ = 0; + for iz = 1:Nz + tmp_ = tmp_ + J(iz)*Napj_(:,:,:,iz,:); + end + Napj_ = tmp_/sum(J(iz)); + % Napj_ = mean(Napj_,4); + % Napj_ = Napj_(:,:,:,Nz/2+1,:); + % phi_ = mean(data.PHI,3); + % otherwise + % iz = options.iz; + % Napj_ = Napj_(:,:,:,:,iz,:); + % phi_ = data.PHI(:,:,iz); +% end % Napj_ = squeeze(Napj_); -frames = options.T; -for it = 1:numel(options.T) - [~,frames(it)] = min(abs(options.T(it)-data.Ts5D)); +frames = T; +for it = 1:numel(T) + [~,frames(it)] = min(abs(T(it)-data.Ts3D)); end frames = unique(frames); -Napj_ = mean(Napj_(:,:,:,:,frames),5); +Napj_ = mean(Napj_(:,:,:,:,frames),5); -% Napj_ = squeeze(Napj_); +Napj_ = squeeze(Napj_); Np = numel(parray); Nj = numel(jarray); -if options.non_adiab - for ij_ = 1:Nj - for ikx = 1:data.Nkx - for iky = 1:data.Nky - kp_ = sqrt(kx_(ikx)^2 + ky_(iky)^2); - Napj_(1,ij_,iky,ikx) = Napj_(1,ij_,iky,ikx) + kernel(ij_,kp_)*phi_(iky,ikx); - end - end - end -end - -if options.RMS - FF = zeros(data.Nky,data.Nkx,numel(options.XPERP),numel(options.SPAR)); +% if options.RMS + FF = zeros(numel(x),numel(s)); FAM = FaM(SS,XX); for ip_ = 1:Np p_ = parray(ip_); HH = Hp(p_,SS); for ij_ = 1:Nj - j_ = jarray(ij_); - LL = Lj(j_,XX); + j_ = jarray(ij_); + LL = Lj(j_,XX); HLF = HH.*LL.*FAM; - for ikx = 1:data.Nkx - for iky = 1:data.Nky - FF(iky,ikx,:,:) = squeeze(FF(iky,ikx,:,:)) + Napj_(ip_,ij_,iky,ikx)*HLF; - end - end + FF = FF + Napj_(ip_,ij_)*HLF; end end -else - FF = zeros(numel(options.XPERP),numel(options.SPAR)); - FAM = FaM(SS,XX); - for ip_ = 1:Np - p_ = parray(ip_); - HH = Hp(p_,SS); - for ij_ = 1:Nj - j_ = jarray(ij_); - LL = Lj(j_,XX); - FF = FF + squeeze(Napj_(ip_,ij_,ikx0,iky0))*HH.*LL.*FAM; - end - end -end +% else +% FF = zeros(numel(options.XPERP),numel(options.SPAR)); +% FAM = FaM(SS,XX); +% for ip_ = 1:Np +% p_ = parray(ip_); +% HH = Hp(p_,SS); +% for ij_ = 1:Nj +% j_ = jarray(ij_); +% LL = Lj(j_,XX); +% FF = FF + squeeze(Napj_(ip_,ij_,ikx0,iky0))*HH.*LL.*FAM; +% end +% end +% end FF = (FF.*conj(FF)); %|f_a|^2 % FF = abs(FF); %|f_a| -if options.RMS +% if options.RMS % FF = squeeze(mean(mean(sqrt(FF),1),2)); %sqrt(<|f_a|^2>kx,ky) - FF = sqrt(squeeze(mean(mean(FF,1),2))); %<|f_a|>kx,ky -else - FF = sqrt(squeeze(FF)); %sqrt(<|f_a|>x,y) -end + FF = sqrt(FF); %<|f_a|>kx,ky +% else +% FF = sqrt(squeeze(FF)); %sqrt(<|f_a|>x,y) +% end -FF = FF./max(max(FF)); -FF = FF'; +% FF = FF./max(max(FF)); +% FF = FF'; % FF = sqrt(FF); % FF = FF'; end \ No newline at end of file diff --git a/matlab/compute/fourier_GENE.m b/matlab/compute/fourier_GENE.m new file mode 100644 index 0000000000000000000000000000000000000000..2039f29ca9ccad0f446dfbe218642b6a97687fac --- /dev/null +++ b/matlab/compute/fourier_GENE.m @@ -0,0 +1,9 @@ +function [ field_c ] = fourier_GENE( field_r ) + +%fft, symmetric as we are using real data +spectrumKxYZ=nx*fft(fieldxyz,[],1); +field_c=ny*fft(spectrumKxYZ,[],2); +clear spectrumKxKyZ + +end + diff --git a/matlab/compute/ifourier_GENE.m b/matlab/compute/ifourier_GENE.m index 1538c36dfb1cde9ae4803a82e3ebf32b7ecddaea..e05abaa4644297f91490770467a25bb20bd35166 100644 --- a/matlab/compute/ifourier_GENE.m +++ b/matlab/compute/ifourier_GENE.m @@ -26,22 +26,6 @@ end spectrumXKyZ=nx*ifft(spectrumKyKxZ,[],1); field_r=ny*ifft(spectrumXKyZ,[],2,'symmetric'); clear spectrumKxKyZ - -%% Adapted for HeLaZ old representation -% [nkx,ny,nz]=size(field_c); -% %extend to whole ky by imposing reality condition -% nx=2*nkx-1; -% -% %note, we need one extra point which we set to zero for the ifft -% spectrumKxKyZ=zeros(nx,ny,nz); -% spectrumKxKyZ(1:nkx,:,:)=field_c(:,:,:); -% spectrumKxKyZ((nkx+1):(nx),1,:)=conj(field_c(nkx:-1:2,1,:)); -% spectrumKxKyZ((nkx+1):(nx),2:ny,:)=conj(field_c(nkx:-1:2,ny:-1:2,:)); -% -% %inverse fft, symmetric as we are using real data -% spectrumXKyZ=nx*ifft(spectrumKxKyZ,[],1); -% field_r=ny*ifft(spectrumXKyZ,[],2,'symmetric'); -% clear spectrumKxKyZ end diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m index ab5736b4f47592164e2158576a322d4c088c1d71..c8998c0a8c22905ed37cf65810a219fa8ccd0fb3 100644 --- a/matlab/compute/mode_growth_meter.m +++ b/matlab/compute/mode_growth_meter.m @@ -6,10 +6,16 @@ d = OPTIONS.fftz.flag; % To add spectral evolution of z (useful for 3d zpin switch OPTIONS.FIELD case 'Ni00' FIELD = DATA.Ni00; - fname = 'N^{00}'; + fname = 'Ni^{00}'; + case 'Ne00' + FIELD = DATA.Ne00; + fname = 'Ne^{00}'; case 'phi' FIELD = DATA.PHI; fname = '\phi'; + case 'T_i' + FIELD = DATA.TEMP_I; + fname = 'T_i'; end if numel(size(FIELD)) == 3 field = squeeze(FIELD); @@ -121,8 +127,13 @@ for i = 1:2 %plot % subplot(2+d,3,1+3*(i-1)) FIGURE.axes(1+3*(i-1)) = subplot(2+d,3,1+3*(i-1),'parent',FIGURE.fig); + if MODES_SELECTOR == 2 [YY,XX] = meshgrid(t,fftshift(k)); pclr = pcolor(XX,YY,abs(plt(fftshift(field,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none'); hold on; + else + [YY,XX] = meshgrid(t,(k)); + pclr = pcolor(XX,YY,abs(plt((field),1:numel(k))));set(pclr, 'edgecolor','none'); hold on; + end set(gca,'YDir','normal') % xlim([t(1) t(end)]); %ylim([1e-5 1]) xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /R_0$'); @@ -138,10 +149,10 @@ for i = 1:2 semilogy(t,plt(field,ikzf),'--k'); end %plot the time window where the gr are measured - plot(t(it1)*[1 1],[1e-10 1],'--k') - plot(t(it2)*[1 1],[1e-10 1],'--k') - xlim([t(1) t(end)]); %ylim([1e-5 1]) + xline(t(it1),'--k') + xline(t(it2),'--k') xlabel('$t c_s /R_0$'); ylabel(['$|',fname,'_{',kstr,'}|$']); + axis tight title('Measure time window') % subplot(2+d,3,3+3*(i-1)) diff --git a/matlab/load/compile_results_2Da.m b/matlab/load/compile_results_2Da.m new file mode 100644 index 0000000000000000000000000000000000000000..f89f413bd464f076b155de47693d5eabb401480b --- /dev/null +++ b/matlab/load/compile_results_2Da.m @@ -0,0 +1,37 @@ +function [field, Ts2D] = compile_results_2Da(DIRECTORY,JOBNUMMIN,JOBNUMMAX,fieldname) +CONTINUE = 1; +JOBNUM = JOBNUMMIN; JOBFOUND = 0; +% field +field = []; +Ts2D = []; +ii = 1; +while(CONTINUE) + filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM); + % Check presence and jobnummax + if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX) + %test if it is corrupted or currently running + try + openable = ~isempty(h5read(filename,'/data/var2d/time')); + catch + openable = 0; + end + if openable + % load field %% + [ F, T, ~] = load_2Da_data(filename, fieldname); + field = cat(4,field,F); + Ts2D = cat(1,Ts2D,T); + ii = ii + 1; + JOBFOUND = JOBFOUND + 1; + end + elseif (JOBNUM > JOBNUMMAX) + CONTINUE = 0; + end + JOBNUM = JOBNUM + 1; +end + +if(JOBFOUND == 0) + disp('no results found, please verify the paths'); + return; +end + +end \ No newline at end of file diff --git a/matlab/load/compile_results_low_mem.m b/matlab/load/compile_results_low_mem.m index 33fd76b5eb6086530d02faf9553cbc2b7dfddc2f..8ef6048ccd5d717c9c2f64699ddbeb84b77c267f 100644 --- a/matlab/load/compile_results_low_mem.m +++ b/matlab/load/compile_results_low_mem.m @@ -29,7 +29,7 @@ while(CONTINUE) DATA.outfilenames{ii} = filename; %test if it is corrupted or currently running try - openable = ~isempty(h5read(filename,'/data/var3d/time')); + openable = ~isempty(h5read(filename,'/data/var0d/time')); catch openable = 0; end @@ -38,13 +38,25 @@ while(CONTINUE) fprintf('Loading ID %.2d (%s)\n',JOBNUM,filename); % loading input parameters DATA.(sprintf('fort_%.2d',JOBNUM)) = struct(); - DATA.(sprintf('fort_%.2d',JOBNUM)) = read_namelist(fortname); + % cp the STDIN.XX from the output file to the directory to load + tmpfort = [fortname,'.tmp']; + fid = fopen(tmpfort,'wt'); + input_string = h5read(filename,sprintf('/files/STDIN.%.2d',JOBNUM)); + fprintf(fid, input_string); + fclose(fid); + DATA.(sprintf('fort_%.2d',JOBNUM)) = read_namelist(tmpfort); + system(['rm ',tmpfort]); % clean it % Loading from output file CPUTIME = h5readatt(filename,'/data/input','cpu_time'); DT_SIM = h5readatt(filename,'/data/input/basic','dt'); [Parray, Jarray, kx, ky, z] = load_grid_data(filename); - W_GAMMA = strcmp(h5readatt(filename,'/data/input/diagnostics','write_gamma'),'y'); - W_HF = strcmp(h5readatt(filename,'/data/input/diagnostics','write_hf' ),'y'); + try + W_GAMMA = strcmp(h5readatt(filename,'/data/input/diagnostics','write_gamma'),'y'); + W_HF = strcmp(h5readatt(filename,'/data/input/diagnostics','write_hf' ),'y'); + catch + W_GAMMA = strcmp(h5readatt(filename,'/data/input/diag_par','write_gamma'),'y'); + W_HF = strcmp(h5readatt(filename,'/data/input/diag_par','write_hf' ),'y'); + end KIN_E = strcmp(h5readatt(filename,'/data/input/model', 'ADIAB_E' ),'n'); BETA = h5readatt(filename,'/data/input/model','beta'); diff --git a/matlab/load/load_2D_data.m b/matlab/load/load_2D_data.m index 67da60c9dd1931e73cde73c7520ecb6370a0c153..6ca82e8430e399ec31c1d9feb74d31e3d2203a6e 100644 --- a/matlab/load/load_2D_data.m +++ b/matlab/load/load_2D_data.m @@ -1,5 +1,5 @@ function [ data, time, dt ] = load_2D_data( filename, variablename ) -%LOAD_3D_DATA load a 3D variable stored in a hdf5 result file from HeLaZ +%LOAD_2D_DATA load a 2D variable stored in a hdf5 result file time = h5read(filename,'/data/var2d/time'); dt = h5readatt(filename,'/data/input/basic','dt'); cstart= h5readatt(filename,'/data/input/basic','start_iframe3d'); diff --git a/matlab/load/load_2Da_data.m b/matlab/load/load_2Da_data.m new file mode 100644 index 0000000000000000000000000000000000000000..3cdd94bfb0f0db936ee9745c89e8755843deea76 --- /dev/null +++ b/matlab/load/load_2Da_data.m @@ -0,0 +1,43 @@ +function [ data, time, dt ] = load_2Da_data( filename, variablename ) +%LOAD_2Da_DATA load a 2D variable stored in a hdf5 result file + time = h5read(filename,'/data/var2d/time'); + dt = h5readatt(filename,'/data/input/basic','dt'); + cstart= h5readatt(filename,'/data/input/basic','start_iframe2d'); + Na = h5readatt(filename,'/data/input/model','Na'); + Np = h5readatt(filename,'/data/input/grid', 'Np'); + Nj = h5readatt(filename,'/data/input/grid', 'Nj'); + Nky = h5readatt(filename,'/data/input/grid', 'Nky'); + Nkx = h5readatt(filename,'/data/input/grid', 'Nkx'); + Nz = h5readatt(filename,'/data/input/grid', 'Nz'); + + + % Find array size by loading the first output + tmp = h5read(filename,['/data/var2d/',variablename,'/', num2str(cstart+1,'%06d')]); + try % check if it is complex or real + sz = size(tmp.real); + cmpx = 1; + catch + sz = size(tmp); + cmpx = 0; + end + % add a dimension if nz=1 or na=1 +% if Na == 1 +% sz = [1 sz]; +% end + if Nz == 1 + sz = [sz 1]; + end + % add time dimension + data = zeros([sz numel(time)]); + + for it = 1:numel(time) + tmp = h5read(filename,['/data/var2d/',variablename,'/', num2str(cstart+it,'%06d')]); + if cmpx + data(:,:,:,it) = reshape(tmp.real,sz) + 1i * reshape(tmp.imaginary,sz); + else + data(:,:,:,it) = reshape(tmp,sz); + end + end + +end + diff --git a/matlab/load/load_multiple_outputs_marconi.m b/matlab/load/load_multiple_outputs_marconi.m deleted file mode 100644 index 1fdd3b8608236ae4a056034abf7852fc50a814a0..0000000000000000000000000000000000000000 --- a/matlab/load/load_multiple_outputs_marconi.m +++ /dev/null @@ -1,6 +0,0 @@ -addpath(genpath('../matlab')) % ... add -load_marconi('/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/') -load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/') -load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/') -load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/500x500_L_120_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_0e+00/') -load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/500x500_L_120_P_2_J_1_eta_0.6_nu_1e-01_DGGK_mu_0e+00/') diff --git a/matlab/load/load_params.m b/matlab/load/load_params.m index 967baf3c4587ab4217c2763ad1c0260becddb13d..85d7c00d0ba3f29c41297abb58c7a6e107cf8b06 100644 --- a/matlab/load/load_params.m +++ b/matlab/load/load_params.m @@ -51,12 +51,19 @@ function [DATA] = load_params(DATA,filename) catch DATA.inputs.ADIAB_I = 0; end - DATA.inputs.W_GAMMA = h5readatt(filename,'/data/input/diagnostics','write_gamma') == 'y'; - DATA.inputs.W_PHI = h5readatt(filename,'/data/input/diagnostics','write_phi') == 'y'; - DATA.inputs.W_NA00 = h5readatt(filename,'/data/input/diagnostics','write_Na00') == 'y'; - DATA.inputs.W_NAPJ = h5readatt(filename,'/data/input/diagnostics','write_Napj') == 'y'; - DATA.inputs.W_SAPJ = h5readatt(filename,'/data/input/diagnostics','write_Sapj') == 'y'; - + try + DATA.inputs.W_GAMMA = h5readatt(filename,'/data/input/diagnostics','write_gamma') == 'y'; + DATA.inputs.W_PHI = h5readatt(filename,'/data/input/diagnostics','write_phi') == 'y'; + DATA.inputs.W_NA00 = h5readatt(filename,'/data/input/diagnostics','write_Na00') == 'y'; + DATA.inputs.W_NAPJ = h5readatt(filename,'/data/input/diagnostics','write_Napj') == 'y'; + DATA.inputs.W_SAPJ = h5readatt(filename,'/data/input/diagnostics','write_Sapj') == 'y'; + catch + DATA.inputs.W_GAMMA = h5readatt(filename,'/data/input/diag_par','write_gamma') == 'y'; + DATA.inputs.W_PHI = h5readatt(filename,'/data/input/diag_par','write_phi') == 'y'; + DATA.inputs.W_NA00 = h5readatt(filename,'/data/input/diag_par','write_Na00') == 'y'; + DATA.inputs.W_NAPJ = h5readatt(filename,'/data/input/diag_par','write_Napj') == 'y'; + DATA.inputs.W_SAPJ = h5readatt(filename,'/data/input/diag_par','write_Sapj') == 'y'; + end % Species dependent parameters DATA.inputs.sigma = zeros(1,DATA.inputs.Na); DATA.inputs.tau = zeros(1,DATA.inputs.Na); diff --git a/matlab/load/readFortranNamelist.m b/matlab/load/readFortranNamelist.m deleted file mode 100644 index 5dbbda9337ba82b2ff2674d9d30573dc51ba30e5..0000000000000000000000000000000000000000 --- a/matlab/load/readFortranNamelist.m +++ /dev/null @@ -1,213 +0,0 @@ -function S = read_namelist(filename) -% S = READ_NAMELIST(FILENAME) returns the struct S containg namelists and -% variables in the file FILENAME organised in hierachical way: -% -% |--VAR1 -% |--VAR2 -% |-- NMLST_A--|... -% | |--VARNa -% | -% | |--VAR1 -% |-- NMLST_B--|--VAR2 -% | |... -% S --| ... |--VARNb -% | -% | |--VAR1 -% |-- NMLST_M--|--VAR2 -% |... -% |--VARNm -% -% Note: The function can read multidimensioal variables as well. The -% function assumes that there is no more than one namelist section per -% line. At this time there is no syntax checking functionality so the -% function will crash in case of errors. -% -% Example: -% NMLST = read_namelist('OPTIONS.nam'); -% NMLST.NAM_FRAC.XUNIF_NATURE = 0.1; -% write_namelist(NMlST, 'MOD_OPTIONS.nam'); -% -% Written by: Darien Pardinas Diaz (darien.pardinas-diaz@monash.edu) -% Version: 1.0 -% Date: 16 Dec 2011 -S = struct(); -% Open and read the text file containing the namelists -fid = fopen(filename,'r'); -c = 0; -lines = cell(1); -% Read all the text lines in namelist file -while ~feof(fid) - line = fgetl(fid); - % Remove comments if any on the line - idx = find(line == '!'); - if ~isempty(idx), - line = line(1:idx(1)-1); - end - if ~isempty(line), - c = c + 1; - lines{c} = line; - end -end -fclose(fid); -i = 0; -while i < c; - % Find a record - i = i + 1; - line = lines{i}; - idx = find(line == '&'); - if ~isempty(idx), % i.e. a namelist start - line = line(idx(1) + 1:end); - % find next space - idx = find(line == ' '); - if ~isempty(idx), - namelst = line(1:idx(1) - 1); - line = line(idx(1) + 1:end); - else - namelst = line; - line = []; - end - nmlst_bdy = []; - idx = strfind(line,'/'); - % Get the variable specification section - while isempty(idx) && i < c, - nmlst_bdy = [nmlst_bdy ' ' line]; - i = i + 1; - line = lines{i}; - idx = strfind(line,'/'); - end - if ~isempty(idx) && idx(1) > 1, - nmlst_bdy = [nmlst_bdy ' ' line]; - end - % Parse current namelist (set of variables) - S.(namelst) = parse_namelist(nmlst_bdy); - end -end -function S = parse_namelist(strng) -% Internal function to parse the body text of a namelist section. -% Limitations: the following patterns are prohibited inside the literal -% strings: '.t.' '.f.' '.true.' '.false.' '(:)' -% Get all .true., .t. and .false., .f. to T and F -strng = regexprep(strng,'\.true\.' ,'T','ignorecase'); -strng = regexprep(strng,'\.false\.','F','ignorecase'); -strng = regexprep(strng,'\.t\.','T','ignorecase'); -strng = regexprep(strng,'\.f\.','F','ignorecase'); -% Make evaluable the (:) expression in MATLAB if any -strng = regexprep(strng, '\(:\)', '(1,:)'); -[strng, islit] = parse_literal_strings([strng ' ']); -% Find the position of all the '=' -eq_idx = find(strng == '='); -nvars = length(eq_idx); -arg_start = eq_idx + 1; -arg_end = zeros(size(eq_idx)); -vars = cell(nvars,1); -S = struct; -% Loop through every variable -for k = 1:nvars, - i = eq_idx(k) - 1; - % Move to the left and discard blank spaces - while strng(i) == ' ', i = i - 1; end - % Now we are over the variable name or closing parentesis - j = i; - if strng(i) == ')', - while strng(i) ~= '(', i = i - 1; end - i = i - 1; - % Move to the left and discard any possible blank spaces - while strng(i) == ' ', i = i - 1; end - end - - % Now we are over the last character of the variable name - while strng(i) ~= ' ', i = i - 1; end - - if k > 1, arg_end(k - 1) = i; end - vars{k} = ['S.' strng(i + 1: j)]; -end -arg_end(end) = length(strng); -% This variables are used in the eval function to evaluate True/False, -% so don't remove it! -T = '.true.'; -F = '.false.'; -% Loop through every variable guess varible type -for k = 1:nvars, - arg = strng(arg_start(k):arg_end(k)); - arglit = islit(arg_start(k):arg_end(k))'; - - % Remove commas in non literal string... - commas = ~arglit & arg == ','; - if any(commas) - arg(commas) = ' '; - end - - if any(arglit), - % We are parsing a variable that is literal string - arg = ['{' arg '};']; - elseif ~isempty(find( arg == 'T' | arg == 'F', 1)), - % We are parsing a boolean variable - arg = ['{' arg '};']; - else - % We are parsing a numerical array - arg = ['[' arg '];']; - end - % Eval the modified syntax in Matlab - eval([vars{k} ' = ' arg]); -end -function [strng, is_lit] = parse_literal_strings(strng) -% Parse the literal declarations of strings and change to Matlab syntax -len = length(strng); -add_squote = []; % Positions to add a scape single quote on syntax -rem_dquote = []; % Positions to remove a double quote scape on syntax -i = 1; -while i < len, - if strng(i) == '''', % Opening string with single quote... - i = i + 1; - while i < len && strng(i) ~= '''' || strcmp(strng(i:i+1),'''''') , - i = i + 1; - if strcmp(strng(i-1:i),''''''), - i = i + 1; - end - end - end - if strng(i) == '"', % Opening string with double quote... - strng(i) = ''''; % Change to single quote - i = i + 1; - while strng(i) ~= '"' || strcmp(strng(i:i+1),'""') && i < len, - % Check for a possible sigle quote here - if strng(i) == '''', - add_squote = [add_squote i]; - end - i = i + 1; - if strcmp(strng(i-1:i),'""'), - rem_dquote = [rem_dquote i-1]; - i = i + 1; - end - end - strng(i) = ''''; % Change to single quote - end - i = i + 1; -end -for i = 1:length(add_squote); - strng = [strng(1:add_squote(i)) strng(add_squote(i):end)]; -end -for i = 1:length(rem_dquote); - strng = [strng(1:rem_dquote(i)-1) strng(rem_squote(i)+1:end)]; -end -% Now everything should be in Matlab string syntax -% Classify syntax as literal or regular expression -i = 1; -len = length(strng); -is_lit = zeros(len,1); -while i < len, - if strng(i) == '''', % Opening string with single quote... - is_lit(i) = 1; - i = i + 1; - while i < len && strng(i) ~= '''' || strcmp(strng(i:i+1),''''''), - is_lit(i) = 1; - i = i + 1; - if strcmp(strng(i-1:i),''''''), - is_lit(i) = 1; - i = i + 1; - end - end - is_lit(i) = 1; - end - i = i + 1; -end diff --git a/matlab/load/read_namelist.m b/matlab/load/read_namelist.m index 3d49e6d0406d2532eaa3fc51ee3f7b2fe84d7631..954dbc9af151a177268b600963ef51968cade8b8 100644 --- a/matlab/load/read_namelist.m +++ b/matlab/load/read_namelist.m @@ -40,10 +40,10 @@ while ~feof(fid) line = fgetl(fid); % Remove comments if any on the line idx = find(line == '!'); - if ~isempty(idx), + if ~isempty(idx) line = line(1:idx(1)-1); end - if ~isempty(line), + if ~isempty(line) c = c + 1; lines{c} = line; end diff --git a/matlab/plot/Ursula_Meier.m b/matlab/plot/Ursula_Meier.m new file mode 100644 index 0000000000000000000000000000000000000000..affabdea29979e8a4f585ad3433beb13caf1df26 --- /dev/null +++ b/matlab/plot/Ursula_Meier.m @@ -0,0 +1,106 @@ +function [filename] = Ursula_Meier(MOVIE) + FPS = 16; BWR = 0; ncolor_level = 128; format = '.gif'; SAT = 0.75; + INTERP = 1; + if isfield(MOVIE,'FPS') + FPS = MOVIE.FPS; + end + if isfield(MOVIE,'BWR') + BWR = MOVIE.BWR; + end + if isfield(MOVIE,'format') + format = MOVIE.format; + end + if isfield(MOVIE,'SAT') + SAT = MOVIE.SAT; + end + if isfield(MOVIE,'INTERP') + INTERP = MOVIE.INTERP; + end + DELAY = 1/MOVIE.FPS; + FILENAME = MOVIE.FILENAME; FIELDNAME = MOVIE.FIELDNAME; + X = MOVIE.X; Y = MOVIE.Y; T = MOVIE.T; F = MOVIE.F; + XNAME = MOVIE.XNAME; YNAME = MOVIE.YNAME; + switch format + case '.avi' + vidfile = VideoWriter(FILENAME,'Uncompressed AVI'); + vidfile.FrameRate = FPS; + open(vidfile); + end + % Set colormap boundaries + hmax = max(max(max(F(:,:,:)))); + hmin = min(min(min(F(:,:,:)))); + if hmax == hmin + disp('Warning : h = hmin = hmax = const') + else + % Setup figure frame + fig = figure('Color','white');%,'Position', toplot.DIMENSIONS.*[0.5 0.5 1.0 1]); + pcolor(X,Y,F(:,:,1)); % to set up + if BWR + colormap(bluewhitered) + else + colormap(gray) + end + clim(SAT*[-1 1]) + axis equal + axis tight % this ensures that getframe() returns a consistent size + if INTERP + shading interp; + end + in = 1; + nbytes = fprintf(2,'frame %d/%d',in,numel(F(1,1,:))); + for n = 1:numel(T) % loop over selected frames + scale = max(max(abs(F(:,:,n)))); % Scaling to normalize + pclr = pcolor(X,Y,F(:,:,n)/scale); % frame plot\ + clim([-1,1]); + if INTERP + shading interp; + end + set(pclr, 'edgecolor','none'); %pbaspect(toplot.ASPECT); + if BWR + colormap(bluewhitered) + else + colormap(gray) + end + clim(SAT*[-1 1]) + title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))... + ,', scaling = ',sprintf('%.1e',scale)]); + + xlabel(XNAME); ylabel(YNAME); %colorbar; + axis tight equal + drawnow + % Capture the plot as an image + frame = getframe(fig); + switch format + case '.gif' + im = frame2im(frame); + [imind,cm] = rgb2ind(im,ncolor_level); + % Write to the GIF File + if in == 1 + imwrite(imind,cm,FILENAME,'gif', 'Loopcount',inf); + else + imwrite(imind,cm,FILENAME,'gif','WriteMode','append', 'DelayTime',DELAY); + end + case '.avi' + writeVideo(vidfile,frame); + otherwise + disp('Unknown format'); + break + end + % terminal info + while nbytes > 0 + fprintf('\b') + nbytes = nbytes - 1; + end + nbytes = fprintf(2,'frame %d/%d',n,numel(F(1,1,:))); + in = in + 1; + end + disp(' ') + switch format + case '.gif' + disp(['Gif saved @ : ',FILENAME]) + case '.avi' + disp(['Video saved @ : ',FILENAME]) + close(vidfile); + end + end +end \ No newline at end of file diff --git a/matlab/plot/cinecita.m b/matlab/plot/cinecita.m new file mode 100644 index 0000000000000000000000000000000000000000..fcfcdd03609bd2a2667102006d974a1b2d372e22 --- /dev/null +++ b/matlab/plot/cinecita.m @@ -0,0 +1,69 @@ +function [ ] = cinecita(DATADIR,J0,J1,fieldname,plan,save) + +data = {}; +[data] = compile_results_low_mem(data,DATADIR,J0,J1); +%% Select the field +SKIP_COMP = 0; % Turned on only for 2D plots +OPE_ = 1; % Operation on data +switch fieldname + case 'phi' %ES pot + [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi'); + ltxname = '\phi'; + case 'phi_obmp' %ES pot + [FIELD,TIME] = compile_results_2D(DATADIR,J0,J1,'phi_obmp'); + ltxname = '\phi_{z=0}'; + SKIP_COMP = 1; + case '\psi' %EM pot + [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'psi'); + ltxname = '\psi'; + case 'phi_nz' % non-zonal ES pot + [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi'); + ltxname = '\phi^{NZ}'; + OPE_ = (KY~=0); + case 'vE_y' % y-comp of ExB velocity + [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi'); + ltxname = 'v_{Ey}'; + OPE_ = -1i*KX; + case 'vE_x' % x-comp of ExB velocity + [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi'); + ltxname = 'v_{Ex}'; + OPE_ = -1i*KY; + case 'sE_y' % y-comp of ExB shear + [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi'); + ltxname = 's_{Ey}'; + OPE_ = KX.^2; + case 'sE_x' % x-comp of ExB shear + [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi'); + ltxname = 's_{Ex}'; + OPE_ = KY.^2; + case 'wz' % ES pot vorticity + [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi'); + ltxname = '\omega_z'; + OPE_ = -(KX.^2+KY.^2); + otherwise + disp('Fieldname not recognized'); + return +end +%% Setup grids +% [KX, KY] = meshgrid(DATA.grids.kx, DATA.grids.ky); +Nx = data.grids.Nx; Ny = data.grids.Ny; Nz = data.grids.Nz; Nt = numel(TIME); +% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Options +options.INTERP = 0; +options.POLARPLOT = 0; +options.BWR = 1; % bluewhitered plot or gray +options.CLIMAUTO = 1; % adjust the colormap auto +options.NAME = ltxname; +options.PLAN = plan; +options.COMP = 1; +options.TIME = data.Time(1:1:end); +data.EPS = 0.1; +data.a = data.EPS * 2000; +options.RESOLUTION = 256; +create_film(data,options,'.gif') + +if ~save + system(['rm ',FILENAME]); +end + +end \ No newline at end of file diff --git a/matlab/plot/cinecita_2D.m b/matlab/plot/cinecita_2D.m new file mode 100644 index 0000000000000000000000000000000000000000..5a6f0f4f5498ac4f3cf3f5a5748ddd1e4857e5e4 --- /dev/null +++ b/matlab/plot/cinecita_2D.m @@ -0,0 +1,136 @@ +function [MOVIE] = cinecita_2D(DATADIR,fieldname,varargin) +DATADIR = [DATADIR,'/']; +J0 = 0; J1 = 20; nskip = 1; FOURIER = 0; SAT = 0.75; +PLAY = 1; +if numel(varargin) > 0 + if isfield(varargin{1},'J0') + J0 = varargin{1}.J0; + end + if isfield(varargin{1},'J1') + J1 = varargin{1}.J1; + end + if isfield(varargin{1},'nskip') + nskip = varargin{1}.nskip; + end + if isfield(varargin{1},'FOURIER') + FOURIER = varargin{1}.FOURIER; + end + if isfield(varargin{1},'SAT') + SAT = varargin{1}.SAT; + end + if isfield(varargin{1},'PLAY') + PLAY = varargin{1}.PLAY; + end +end +% +data = {}; +data = compile_results_low_mem(data,DATADIR,J0,J1); +data = load_params(data,data.outfilenames{1}); +Nx = data.grids.Nx; Ny = data.grids.Ny; Na = data.inputs.Na; +[KX,KY] = meshgrid(data.grids.kx,data.grids.ky); +format = '.gif'; +%% Select the field +OPE_ = 1; % Operation on data +Nspec = 0; % Load a specie +switch fieldname +case '\phi' %ES pot + toload = 'phi'; + ltxname = '\phi'; + vname = 'phi'; +case '\phi_{zf}' %ES pot + toload = 'phi'; + OPE_ = KY==0; % Operation on data + ltxname = '\phi_{zf}'; + vname = 'phizf'; +case '\psi' %EM pot + toload = 'psi'; + ltxname = '\psi'; + vname = 'psi'; +case 'n_i' % ion density + toload = 'dens'; + Nspec = 1; + ltxname = 'n_i'; + vname = 'dens_i'; +case 'n_e' % ion density + toload = 'dens'; + Nspec = 2; + ltxname = 'n_e'; + vname = 'dens_e'; +case 'v_{Ex}' %ES pot + toload = 'phi'; + OPE_ = -1i*KY; % Operation on data + ltxname= 'v_{Ex}'; + vname = 'vEx'; +case 'v_{Ey}' %ES pot + toload = 'phi'; + OPE_ = 1i*KX; % Operation on data + ltxname= 'v_{Ey}'; + vname = 'vEy'; +case '\delta B_x' + toload = 'psi'; + OPE_ = -1i*KY; % Operation on data + ltxname = '\delta B_x'; + vname = 'dBx'; +case '\delta B_y' + toload = 'psi'; + OPE_ = 1i*KX; % Operation on data + ltxname = '\delta B_y'; + vname = 'dBy'; +otherwise + toload = fieldname; + OPE_ = 1; +end + +if Nspec + [FIELD,TIME] = compile_results_2Da(DATADIR,J0,J1,[toload,'_obmp']); + FIELD = squeeze(FIELD(Nspec,:,:,:)); +else + [FIELD,TIME] = compile_results_2D(DATADIR,J0,J1,[toload,'_obmp']); +end + +TIME = TIME(1:nskip:end); FIELD = FIELD(:,:,1:nskip:end); + +switch ~FOURIER + case 1 % Real space plot + XNAME = '$x$'; YNAME = '$y$'; plane ='xy'; + [X,Y] = meshgrid(data.grids.x,data.grids.y); + INTERP = 1; + process = @(x) real(fftshift(ifourier_GENE(x))); + shift_x = @(x) x; + shift_y = @(x) x; + case 0 % Frequencies plot + XNAME = '$k_x$'; YNAME = '$k_y$'; plane ='kxky'; + [X,Y] = meshgrid(data.grids.kx,data.grids.ky); + INTERP = 0; + process = @(x) abs(fftshift(x,2)); + shift_x = @(x) fftshift(x,2); + shift_y = @(x) fftshift(x,2); +end + +for it = 1:numel(TIME) + tmp = squeeze(OPE_.*FIELD(:,:,it)); + F(:,:,it) = squeeze(process(tmp)); +end + +X = shift_x(X); +Y = shift_y(Y); +FILENAME = [vname,'_',plane,'_obmp']; +FILENAME = [DATADIR,FILENAME,format]; +FIELDNAME = ['$',ltxname,'$']; +MOVIE.F = F; +MOVIE.X = X; +MOVIE.Y = Y; +MOVIE.T = TIME; +MOVIE.XNAME = XNAME; +MOVIE.YNAME = YNAME; +MOVIE.FPS = 16; +MOVIE.BWR = 0; +MOVIE.SAT = 0.75; +MOVIE.FIELDNAME = FIELDNAME; +MOVIE.FILENAME = FILENAME; +% Shoot the movie +if PLAY + Ursula_Meier(MOVIE) +end + +end \ No newline at end of file diff --git a/matlab/plot/create_film.m b/matlab/plot/create_film.m index 6becccec8c45ab7c06a95c19e033bbf3d28363c9..ca56570d89752bf52d946ac16a07bc43a0a40f99 100644 --- a/matlab/plot/create_film.m +++ b/matlab/plot/create_film.m @@ -1,16 +1,22 @@ function create_film(DATA,OPTIONS,format) %% Plot options -FPS = 16; DELAY = 1/FPS; +FPS = OPTIONS.FPS; DELAY = 1/FPS; BWR = OPTIONS.BWR; NORMALIZED = 1; -if ~strcmp(OPTIONS.PLAN,'sx') - T = DATA.Ts3D; -else - T = DATA.Ts5D; -end %% Processing switch OPTIONS.PLAN case 'RZ' toplot = poloidal_plot(DATA,OPTIONS); + case '3D' + OPTIONS.PLAN = 'xy'; + [~,OPTIONS.COMP] = min(abs(OPTIONS.XYZ(3) - DATA.grids.z)); + toplot = process_field(DATA,OPTIONS); + OPTIONS.PLAN = 'xz'; + [~,OPTIONS.COMP] = min(abs(OPTIONS.XYZ(2) - DATA.grids.y)); + toplot_xz = process_field(DATA,OPTIONS); + OPTIONS.PLAN = 'yz'; + [~,OPTIONS.COMP] = min(abs(OPTIONS.XYZ(1) - DATA.grids.x)); + toplot_yz = process_field(DATA,OPTIONS); + OPTIONS.PLAN = '3D'; otherwise toplot = process_field(DATA,OPTIONS); end @@ -20,7 +26,7 @@ XNAME = ['$',toplot.XNAME,'$']; YNAME = ['$',toplot.YNAME,'$']; FIELDNAME = ['$',toplot.FIELDNAME,'$']; FIELD = toplot.FIELD; X = toplot.X; Y = toplot.Y; -FRAMES = toplot.FRAMES; +TIME = toplot.TIME; switch format case '.avi' vidfile = VideoWriter(FILENAME,'Uncompressed AVI'); @@ -36,9 +42,22 @@ if hmax == hmin disp('Warning : h = hmin = hmax = const') else % Setup figure frame -fig = figure('Color','white');%,'Position', toplot.DIMENSIONS.*[0.5 0.5 1.0 1]); if ~strcmp(OPTIONS.PLAN,'sx') - pcolor(X,Y,FIELD(:,:,1)); % to set up + if ~strcmp(OPTIONS.PLAN,'3D') + fig = figure('Color','white'); + pclr = pcolor(X,Y,FIELD(:,:,1)/scale); % frame plot\ + set(pclr, 'edgecolor','none'); %pbaspect(toplot.ASPECT); + else + fig = figure('Color','white','Position', 1024*[0.5 0.5 0.5 1]); + s = surface(toplot.X,toplot.Y,OPTIONS.XYZ(3)+0*toplot.X,FIELD(:,:,1)/scale); hold on; + s.EdgeColor = 'none'; + s = surface(toplot_xz.X,OPTIONS.XYZ(2)+0*toplot_xz.X,toplot_xz.Y,toplot_xz.FIELD(:,:,1)./scale); + s.EdgeColor = 'none'; + s = surface(OPTIONS.XYZ(1)+0*toplot_yz.X,toplot_yz.X,toplot_yz.Y,toplot_yz.FIELD(:,:,1)./scale); + s.EdgeColor = 'none'; + zlabel('z'); + view([1 -1 0.25]) + end if BWR colormap(bluewhitered) else @@ -52,30 +71,31 @@ fig = figure('Color','white');%,'Position', toplot.DIMENSIONS.*[0.5 0.5 1.0 1]) shading interp; end else + fig = figure('Color','white'); contour(toplot.X,toplot.Y,FIELD(:,:,1),128); end in = 1; nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:))); - for n = 1:numel(FRAMES) % loop over selected frames + for n = 1:numel(TIME) % loop over selected frames scale = max(max(abs(FIELD(:,:,n)))); % Scaling to normalize if ~strcmp(OPTIONS.PLAN,'sx') - if NORMALIZED + if ~strcmp(OPTIONS.PLAN,'3D') pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot\ - caxis([-1,1]); + set(pclr, 'edgecolor','none'); %pbaspect(toplot.ASPECT); else - pclr = pcolor(X,Y,FIELD(:,:,n)); % frame plot\ - if CONST_CMAP - caxis([-1,1]*max(abs([hmin hmax]))); % adaptive color map - else - caxis([-1,1]*scale); % adaptive color map - end - title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))... - ,', scale = ',sprintf('%.1e',scale)]); + s = surface(toplot.X,toplot.Y,OPTIONS.XYZ(3)+0*toplot.X,FIELD(:,:,n)/scale); hold on; + s.EdgeColor = 'none'; + s = surface(toplot_xz.X,OPTIONS.XYZ(2)+0*toplot_xz.X,toplot_xz.Y,toplot_xz.FIELD(:,:,n)./scale); + s.EdgeColor = 'none'; + s = surface(OPTIONS.XYZ(1)+0*toplot_yz.X,toplot_yz.X,toplot_yz.Y,toplot_yz.FIELD(:,:,n)./scale); + s.EdgeColor = 'none'; + zlabel('z'); + view([1 -1 0.25]) end + clim([-1,1]); if toplot.INTERP shading interp; end - set(pclr, 'edgecolor','none'); %pbaspect(toplot.ASPECT); if BWR colormap(bluewhitered) else @@ -87,10 +107,14 @@ fig = figure('Color','white');%,'Position', toplot.DIMENSIONS.*[0.5 0.5 1.0 1]) else % show velocity distr. contour(toplot.X,toplot.Y,FIELD(:,:,n)/scale,128); end - title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(FRAMES(n))))... - ,', scaling = ',sprintf('%.1e',scale)]); - - xlabel(XNAME); ylabel(YNAME); %colorbar; + if ~OPTIONS.RMAXIS + title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(TIME(n)))... + ,', scaling = ',sprintf('%.1e',scale)]); + xlabel(XNAME); ylabel(YNAME); %colorbar; + else + axis off + axis tight + end drawnow % Capture the plot as an image frame = getframe(fig); diff --git a/matlab/plot/photomaton.m b/matlab/plot/photomaton.m index 67b02747b9ca4d0f12ce50d12ff2a476018ecc3e..d563030d9b1ae66ad1ab1cb54ac69dec98bd3b4d 100644 --- a/matlab/plot/photomaton.m +++ b/matlab/plot/photomaton.m @@ -4,11 +4,26 @@ function [ FIGURE ] = photomaton( DATA,OPTIONS ) switch OPTIONS.PLAN case 'RZ' toplot = poloidal_plot(DATA,OPTIONS); + case '3D' + OPTIONS.PLAN = 'xy'; + [~,OPTIONS.COMP] = min(abs(OPTIONS.XYZ(3) - DATA.grids.z)); + toplot = process_field(DATA,OPTIONS); + OPTIONS.PLAN = 'xz'; + [~,OPTIONS.COMP] = min(abs(OPTIONS.XYZ(2) - DATA.grids.y)); + toplot_xz = process_field(DATA,OPTIONS); + OPTIONS.PLAN = 'yz'; + [~,OPTIONS.COMP] = min(abs(OPTIONS.XYZ(1) - DATA.grids.x)); + toplot_yz = process_field(DATA,OPTIONS); + OPTIONS.PLAN = '3D'; otherwise toplot = process_field(DATA,OPTIONS); end FNAME = toplot.FILENAME; FRAMES = toplot.FRAMES; +if OPTIONS.TAVG + toplot.FIELD = mean(toplot.FIELD,3); + FRAMES = FRAMES(1); +end Nframes= numel(FRAMES); Nrows = ceil(Nframes/4); Ncols = ceil(Nframes/Nrows); @@ -29,17 +44,31 @@ FIGURE.fig = figure; %set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows] end if ~strcmp(OPTIONS.PLAN,'sx') tshot = DATA.Ts3D(FRAMES(i_)); - pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale); set(pclr, 'edgecolor','none'); - if OPTIONS.AXISEQUAL - pbaspect(toplot.ASPECT) - end - if ~strcmp(OPTIONS.PLAN,'kxky') - clim([-1,1]*frame_max/scale); - colormap(bluewhitered); - if OPTIONS.INTERP - shading interp; - end - end + if ~strcmp(OPTIONS.PLAN,'3D') + pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale); set(pclr, 'edgecolor','none'); + else + % xy plane + s = surface(toplot.X,toplot.Y,OPTIONS.XYZ(3)+0*toplot.Y,toplot.FIELD(:,:,i_)./scale); hold on; + s.EdgeColor = 'none'; + % s.FaceAlpha = 0.5; + % xz plane + s = surface(toplot_xz.X,OPTIONS.XYZ(2)+0*toplot_xz.X,toplot_xz.Y,toplot_xz.FIELD(:,:,i_)./scale); + s.EdgeColor = 'none'; + % s.FaceAlpha = 0.7; + % yz plane + s = surface(OPTIONS.XYZ(1)+0*toplot_yz.X,toplot_yz.X,toplot_yz.Y,toplot_yz.FIELD(:,:,i_)./scale); + s.EdgeColor = 'none'; + % s.FaceAlpha = 0.7; + xlabel('x'); + ylabel('y'); + zlabel('z'); + h = gca; + plot3(h.XLim,OPTIONS.XYZ(2)*[1 1],OPTIONS.XYZ(3)*[1 1],'--k','LineWidth',1.) + plot3(OPTIONS.XYZ(1)*[1 1],h.YLim,OPTIONS.XYZ(3)*[1 1],'--k','LineWidth',1.) + plot3(OPTIONS.XYZ(1)*[1 1],OPTIONS.XYZ(2)*[1 1],h.ZLim,'--k','LineWidth',1.) + % zline(OPTIONS.XYZ(3),'-k','LineWidth',1.5) + view([1 -1 0.25]) + end else contour(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale,128); % pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale); set(pclr, 'edgecolor','none'); shading interp @@ -47,12 +76,23 @@ FIGURE.fig = figure; %set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows] end xlabel(toplot.XNAME); ylabel(toplot.YNAME); % if i_ > 1; set(gca,'ytick',[]); end; - title([sprintf('$t c_s/R=%.0f$',tshot),', max = ',sprintf('%.1e',frame_max)]); + title([sprintf('$t c_s/R=%5.2f$',tshot),', max = ',sprintf('%.1e',frame_max)]); if OPTIONS.CLIMAUTO clim('auto') end end + if OPTIONS.AXISEQUAL + pbaspect(toplot.ASPECT) + end + if ~strcmp(OPTIONS.PLAN,'kxky') + clim([-1,1]*frame_max/scale); + colormap(bluewhitered); + if OPTIONS.INTERP + shading interp; + end + end legend(['$',toplot.FIELDNAME,'$']); + legend('Location','northeast'); FIGURE.FIGNAME = [FNAME,'_snaps',TNAME]; end diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m index f43866e82dd9cc6cf5cb02dad988bbd859b9f083..aa6565c71b034dbd4be27b87d4be6e42c20b534d 100644 --- a/matlab/plot/plot_ballooning.m +++ b/matlab/plot/plot_ballooning.m @@ -6,12 +6,12 @@ function [FIG, b_angle, phib, psib, ky_] = plot_ballooning(DATA,OPTIONS) [~,ikyarray] = min(abs(DATA.grids.ky - OPTIONS.kymodes)); ikyarray = unique(ikyarray); dz = DATA.grids.z(2) - DATA.grids.z(1); - phi_real=real(DATA.PHI(:,:,:,it1)); - phi_imag=imag(DATA.PHI(:,:,:,it1)); + phi_real=mean(real(DATA.PHI(:,:,:,it0:it1)),4); + phi_imag=mean(imag(DATA.PHI(:,:,:,it0:it1)),4); ncol = 1; if DATA.inputs.BETA > 0 - psi_real=real(DATA.PSI(:,:,:,it1)); - psi_imag=imag(DATA.PSI(:,:,:,it1)); + psi_real=mean(real(DATA.PSI(:,:,:,it0:it1)),4); + psi_imag=mean(imag(DATA.PSI(:,:,:,it0:it1)),4); ncol = 2; end if OPTIONS.PLOT_KP @@ -50,12 +50,12 @@ function [FIG, b_angle, phib, psib, ky_] = plot_ballooning(DATA,OPTIONS) title(['$k_y=',sprintf('%2.2f',DATA.grids.ky(iky)),... ',t_{avg}\in [',sprintf('%1.1f',DATA.Ts3D(it0)),',',sprintf('%1.1f',DATA.Ts3D(it1)),']$']); % z domain start end point - for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1) - xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on; - end - for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1) - xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on; - end + % for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1) + % xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on; + % end + % for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1) + % xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on; + % end if DATA.inputs.BETA > 0 [psib_real,b_angle] = ballooning_representation(psi_real,DATA.inputs.SHEAR,DATA.grids.Npol,DATA.grids.kx,iky,DATA.grids.ky,DATA.grids.z); [psib_imag,~] = ballooning_representation(psi_imag,DATA.inputs.SHEAR,DATA.grids.Npol,DATA.grids.kx,iky,DATA.grids.ky,DATA.grids.z); @@ -71,12 +71,12 @@ function [FIG, b_angle, phib, psib, ky_] = plot_ballooning(DATA,OPTIONS) legend('real','imag','norm') % z domain start end point - for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1) - xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on; - end - for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1) - xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on; - end + % for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1) + % xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on; + % end + % for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1) + % xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on; + % end xlabel('$\chi / \pi$') ylabel('$\psi_B (\chi)$'); title(['$k_y=',sprintf('%2.2f',DATA.grids.ky(iky)),... @@ -89,12 +89,12 @@ function [FIG, b_angle, phib, psib, ky_] = plot_ballooning(DATA,OPTIONS) subplot(numel(ikyarray),ncol,ncol*(iplot-1)+3) plot(b_angle/pi, kperpb,'-k'); hold on; % z domain start end point - for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1) - xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on; - end - for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1) - xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on; - end + % for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1) + % xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on; + % end + % for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1) + % xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on; + % end xlabel('$\chi / \pi$') ylabel('$k_\perp (\chi)$'); title(['$k_y=',sprintf('%2.2f',DATA.grids.ky(iky)),... diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m index ff31e4d26519c7a4f4f2870216fdd478d01a2698..3954ced9b260fe2f729e172b68e32cddb4c0a6c4 100644 --- a/matlab/plot/plot_cosol_mat.m +++ b/matlab/plot/plot_cosol_mat.m @@ -1,97 +1,103 @@ -% 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'); +%% Coeff trajectory +% matfilename = 'gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5'; JMAX = 5; +% matfilename = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'; JMAX = 10; +% matfilename = 'gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'; JMAX = 10; +matfilename = 'gk_landau_P16_J9_dk_5e-2_km_2.0_NFLR_8.h5'; JMAX = 10; +P_ = 10; J_ = P_/2; +n = 1:((P_+1)*(J_+1)); +for p_ = 0:P_ + n((0:J_)+p_*(J_+1)+1) = (0:J_)+p_*(JMAX+1); +end +mat_file_name = ['/home/ahoffman/gyacomo/iCa/',matfilename]; +kp_a = h5read(mat_file_name,'/coordkperp'); +coeff = kp_a*0; +gmax = kp_a*0; +p1 = 0; j1 = 0; +p2 = 0; j2 = 0; +for ik = 1:numel(kp_a) + matidx = sprintf('%5.5i',ik-1); + M_self = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']); + MAT = M_self(n+1,n+1); + ipj1 = n(j1+p1*(J_+1)+1); + ipj2 = n(j2+p2*(J_+1)+1); + coeff(ik) = MAT(ipj1+1,ipj2+1); + gmax(ik) = -max((real(eig(MAT)))); + +end +figure +plot(kp_a,gmax); + %% SGGK +P_ = 10; J_ = 5; +n = 1:((P_+1)*(J_+1)); +for p_ = 0:P_ + n((0:J_)+p_*(J_+1)+1) = (0:J_)+p_*(10+1); +end +kp = 5.0; +kp_a = h5read(mat_file_name,'/coordkperp'); +[~,matidx] = min(abs(kp_a-kp)); +matidx = sprintf('%5.5i',matidx); +mat_file_name = '/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'; +M_self = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']); 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'); +MAT = M_self; %suptitle('SGGK ii,P=20,J=10, k=0'); +imagesc((MAT(n+1,n+1))); +title('Sugama') +xlim('tight') +ylim('tight') +axis equal +% clim(1*[-1 1]) +clim auto +colormap(bluewhitered) %% 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'); - +P_ = 10; J_ = 5; +n = 1:((P_+1)*(J_+1)); +for p_ = 0:P_ + n((0:J_)+p_*(J_+1)+1) = (0:J_)+p_*(10+1); +end +kp = 1.0; +kp_a = h5read(mat_file_name,'/coordkperp'); +[~,matidx] = min(abs(kp_a-kp)); +matidx = sprintf('%5.5i',matidx); +mat_file_name = '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'; +PAGK_self = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']); 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'); +PA_MAT = PAGK_self; + imagesc((PA_MAT(n+1,n+1))); + title('Lorentz') + xlim('tight') + ylim('tight') + axis equal + % clim(1*[-1 1]) + clim auto + colormap(bluewhitered) %% FCGK -P_ = 4; J_ = 2; -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 = 0.0; +P_ = 10; J_ = 5; +n = 1:((P_+1)*(J_+1)); +for p_ = 0:P_ + n((0:J_)+p_*(J_+1)+1) = (0:J_)+p_*(5+1); +end +mat_file_name = '/home/ahoffman/gyacomo/iCa/gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5'; +kp = 1.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,')']); +FC_MAT = FCGK_self; % 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'); + imagesc((FC_MAT(n+1,n+1))); + title('Landau') + xlim('tight') + ylim('tight') + axis equal + % clim(1*[-1 1]) + clim auto + colormap(bluewhitered) %% Eigenvalue spectrum analysis if 0 @@ -129,16 +135,19 @@ for j_ = 1:numel(mfns) % kp_a = kp_a(kp_a<=3); gmax = zeros(size(kp_a)); wmax = zeros(size(kp_a)); + c44 = 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((real(eig(MAT)))); - wmax(idx_+1) = max((imag(eig(MAT)))); + FC_MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']); + gmax(idx_+1) = max((real(eig(FC_MAT)))); + wmax(idx_+1) = max((imag(eig(FC_MAT)))); + c44 (idx_+1) = FC_MAT(4,4); 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; + % plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on; + plot(kp_a,c44,'DisplayName',CONAME_A{j_}); hold on; end subplot(121) legend('show'); grid on; @@ -195,9 +204,9 @@ for j_ = 1:numel(mfns) 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))*TAU_A(j_); - puls{j_} = imag(eig(MAT))*TAU_A(j_); + FC_MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']); + grow{j_} = real(eig(FC_MAT))*TAU_A(j_); + puls{j_} = imag(eig(FC_MAT))*TAU_A(j_); end figure diff --git a/matlab/plot/plot_isosurf_wip.m b/matlab/plot/plot_isosurf_wip.m new file mode 100644 index 0000000000000000000000000000000000000000..886c29fa643311a250986760a1197105ab3b7de7 --- /dev/null +++ b/matlab/plot/plot_isosurf_wip.m @@ -0,0 +1,20 @@ +options.NAME = ['N_i^{00}']; +options.PLAN = 'xyz'; +options.TIME = [30]; +TOPLOT = process_field( data, options ); + +[X,Y,Z] = meshgrid(data.grids.x,data.grids.y,data.grids.z); +X = smooth3(X,'gaussian',5); +Y = smooth3(Y,'gaussian',5); +Z = smooth3(Z,'gaussian',5); +TOPLOT.FIELD = smooth3(TOPLOT.FIELD,'gaussian',5); +s = isosurface(X,Y,Z,TOPLOT.FIELD,0); +s = reducepatch(s,0.1); +%% figure + +figure + +p=patch(s); +set(p,'EdgeColor','none'); +set(p,'FaceColor','r'); +set(p,'FaceAlpha',0.4); diff --git a/matlab/plot/plot_kernels.m b/matlab/plot/plot_kernels.m index a02b82cae6367eb22d4012c19859286eec4407cb..2e7f700e9878c13d7766329bfc8f7bf4f915d86b 100644 --- a/matlab/plot/plot_kernels.m +++ b/matlab/plot/plot_kernels.m @@ -2,15 +2,15 @@ if 0 %% Kernels kmax=5; nmax=6; -kx = linspace(0,kmax,100); +kp = linspace(0,kmax,100); figure for n_ = 0:nmax - plot(kx,kernel(n_,kx),'DisplayName',['$\mathcal{K}_{',num2str(n_),'}$']);hold on; + plot(kp,kernel(n_,kp),'DisplayName',['$\mathcal{K}_{',num2str(n_),'}$']);hold on; end ylim_ = ylim; -plot(kx(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA'); +plot(kp(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA'); % J0 = besselj(0,kx); % plot(kx,J0,'-r','DisplayName','$J_0$'); legend('show'); xlabel('k'); title('Kernel functions $\mathcal{K}_n$'); @@ -19,32 +19,32 @@ end if 0 %% Bessels and approx -vperp = linspace(0,1.5,4); +wperp = linspace(0,1.5,4); nmax1=5; nmax2=10; kmax=7; figure -for i = 1:4 -subplot(2,2,i) - v_ = vperp(i); - kx = linspace(0,kmax,100); +for id2 = 1:4 +subplot(2,2,id2) + v_ = wperp(id2); + kp = linspace(0,kmax,100); - J0 = besselj(0,kx*v_); - A1 = 1 - kx.^2*v_^2/4; - K1 = zeros(size(kx)); - K2 = zeros(size(kx)); + J0 = besselj(0,kp*v_); + A1 = 1 - kp.^2*v_^2/4; + K1 = zeros(size(kp)); + K2 = zeros(size(kp)); for n_ = 0:nmax1 - K1 = K1 + kernel(n_,kx).*polyval(LaguerrePoly(n_),v_^2); + K1 = K1 + kernel(n_,kp).*polyval(LaguerrePoly(n_),v_^2); end for n_ = 0:nmax2 - K2 = K2 + kernel(n_,kx).*polyval(LaguerrePoly(n_),v_^2); + K2 = K2 + kernel(n_,kp).*polyval(LaguerrePoly(n_),v_^2); end - plot(kx,J0,'-k','DisplayName','$J_0(kv)$'); hold on; - plot(kx,A1,'-r','DisplayName','$1 - k^2 v^2/4$'); - plot(kx,K1,'--b','DisplayName',['$\sum_{n=0}^{',num2str(nmax1),'}\mathcal K_n(k) L^n(v)$']); - plot(kx,K2,'-b','DisplayName',['$\sum_{n=0}^{',num2str(nmax2),'}\mathcal K_n(k) L^n(v)$']); + plot(kp,J0,'-k','DisplayName','$J_0(kv)$'); hold on; + plot(kp,A1,'-r','DisplayName','$1 - k^2 v^2/4$'); + plot(kp,K1,'--b','DisplayName',['$\sum_{n=0}^{',num2str(nmax1),'}\mathcal K_n(k) L^n(v)$']); + plot(kp,K2,'-b','DisplayName',['$\sum_{n=0}^{',num2str(nmax2),'}\mathcal K_n(k) L^n(v)$']); ylim_ = [-0.5, 1.0]; - plot(kx(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA'); + plot(kp(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA'); ylim(ylim_); xlabel('$k$') legend('show'); grid on; title(['$v = ',num2str(v_),'$']) end @@ -53,165 +53,223 @@ end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if 0 %% from Sum_n Kn x Ln x Lj to Sum_n Kn Sum_s dnjs x Ls -vperp = [2]; -kx = linspace(0,10,200); -Jmax = 15; -j = 10; -fig = figure; set(gcf, 'Position', [100, 100, numel(vperp)*300, 250]) -% suptitle(['$J_{max}=',num2str(Jmax),', j=',num2str(j),'$']) +wperp = [0 0.2 0.4]; +kp = linspace(0,3,50); +Jmax = 8; +js = 4; +fig = figure; set(gcf, 'Position', [100, 100, numel(wperp)*300, 250]) +% suptitle(['$J_{ma}=',num2str(Jmax),', j=',num2str(j),'$']) Kn = @(x__,y__) kernel(x__,y__); -for i = 1:numel(vperp) -subplot(1,numel(vperp),i) - v_ = vperp(i); +% dnjs_ = @(n,j,s) dnjs(n,j,s); +dnjs_ = @(n,j,s) dnjs_GYAC(n+1,j+1,s+1); +for id2 = 1:numel(wperp) +subplot(1,numel(wperp),id2) + v_ = wperp(id2); % J0 x Lj - J0Lj = besselj(0,kx*v_)*polyval(LaguerrePoly(j),v_^2); + J0Lj_ = besselj(0,kp*v_)*polyval(LaguerrePoly(js),v_^2); % Taylor exp of J0 - A1 = (1 - kx.^2*v_^2/4)*polyval(LaguerrePoly(j),v_^2); + A1 = (1 - kp.^2*v_^2/4)*polyval(LaguerrePoly(js),v_^2); % Sum_n^Nmax Kn x Ln x Lj - KLL = zeros(size(kx)); + KLL = zeros(size(kp)); for n_ = 0:Jmax - KLL = KLL + Kn(n_,kx).*polyval(LaguerrePoly(n_),v_^2); + KLL = KLL + Kn(n_,kp).*polyval(LaguerrePoly(n_),v_^2); end - KLL = KLL.*polyval(LaguerrePoly(j),v_^2); + KLL = KLL.*polyval(LaguerrePoly(js),v_^2); % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls - KdL1 = zeros(size(kx)); - for n_ = 0:Jmax-j + KdL1_ = zeros(size(kp)); + for n_ = 0:Jmax-js tmp_ = 0; - for s_ = 0:n_+j - tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2); + for s_ = 0:n_+js + tmp_ = tmp_ + dnjs_(n_,js,s_)*polyval(LaguerrePoly(s_),v_^2); end - KdL1 = KdL1 + Kn(n_,kx).*tmp_; + KdL1_ = KdL1_ + Kn(n_,kp).*tmp_; end % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls - KdL2 = zeros(size(kx)); + KdL2_ = zeros(size(kp)); for n_ = 0:Jmax tmp_ = 0; - for s_ = 0:min(Jmax,n_+j) - tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2); + for s_ = 0:min(Jmax,n_+js) + tmp_ = tmp_ + dnjs_(n_,js,s_)*polyval(LaguerrePoly(s_),v_^2); end - KdL2 = KdL2 + Kn(n_,kx).*tmp_; + KdL2_ = KdL2_ + Kn(n_,kp).*tmp_; end % plots - plot(kx,J0Lj,'-k','DisplayName','$J_0 L_j$'); hold on; - plot(kx,KLL,'-','DisplayName',['$\sum_{n=0}^{J}\mathcal K_n L^n L_j$']); - plot(kx,KdL1,'-.','DisplayName',['$\sum_{n=0}^{J-j}\mathcal K_n \sum_{s=0}^{n+j} d_{njs} L^s$']); + plot(kp,J0Lj_,'-k','DisplayName','$J_0 L_j$'); ylim_ = ylim; - plot(kx,A1,'-','DisplayName','$(1 - k^2 v^2/4) L_j$'); hold on; - plot(kx,KdL2,'--','DisplayName',['$\sum_{n=0}^{J}\mathcal K_n \sum_{s=0}^{\min(J,n+j)} d_{njs} L^s$']); + plot(kp,KLL,'-','DisplayName',['$\sum_{n=0}^{J}\mathcal K_n L_n L_j$']); hold on; + plot(kp,KdL2_,':o','DisplayName',['$\sum_{n=0}^{J}\mathcal K_n \sum_{s=0}^{\min(J,n+j)} d_{njs} L_s$']); + plot(kp,J0Lj_,'-k','DisplayName','$J_0 L_j$'); hold on; + ylim_ = ylim; + plot(kp,A1,':s','DisplayName','$(1 - l_\perp) L_j$'); hold on; + plot(kp,KdL1_,':x','MarkerSize',10,... + 'DisplayName',['$\sum_{n=0}^{J-j}\mathcal K_n \sum_{s=0}^{n+j} d_{njs} L_s$']); + plot(kp,J0Lj_,'-k','DisplayName','$J_0 L_j$'); hold on; % plot(kx(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA'); ylim(ylim_); - xlabel('$k_{\perp}$') - legend('show'); - grid on; title(['$v = ',num2str(v_),'$',' $J_{max}=',num2str(Jmax),', j=',num2str(j),'$']) + xlabel('$k_{\perp}\rho_s$') + % legend('show'); + grid on; title(['$w_\perp = ',num2str(v_),'$',' $J_{max}=',num2str(Jmax),', j=',num2str(js),'$']) end -FIGNAME = ['/home/ahoffman/Dropbox/Applications/Overleaf/Hermite-Laguerre Z-pinch (HeLaZ Doc.)/figures/J0Lj_approx_J','_',num2str(Jmax),'_j_',num2str(j),'.eps']; +FIGNAME = ['/home/ahoffman/Dropbox/Applications/Overleaf/Hermite-Laguerre Z-pinch (HeLaZ Doc.)/figures/J0Lj_approx_J','_',num2str(Jmax),'_j_',num2str(js),'.eps']; % saveas(fig,FIGNAME,'epsc'); end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -if 1 -%% Convergence analysis -kx = linspace(0,10,200); -v_A = linspace(0,1,10); -J_A = 1:15; -[YY,XX] = meshgrid(v_A,J_A); -klimLL = zeros(size(XX)); -klimdL1= zeros(size(XX)); -Kn = @(x__,y__) kernel(x__,y__); -fig = figure; set(gcf, 'Position', [100, 100, 500, 400]); -for j = 5 -for ij_ = 1:numel(J_A) - iv_ = 1; - for v_ = v_A - eps = abs(0.01*polyval(LaguerrePoly(j),v_^2)); - Jmax = J_A(ij_); - % J0 x Lj - J0Lj = besselj(0,kx*v_)*polyval(LaguerrePoly(j),v_^2); - % Taylor exp of J0 - A1 = (1 - kx.^2*v_^2/4)*polyval(LaguerrePoly(j),v_^2); - % Sum_n^Nmax Kn x Ln x Lj - KLL = zeros(size(kx)); - for n_ = 0:Jmax - KLL = KLL + Kn(n_,kx).*polyval(LaguerrePoly(n_),v_^2); - end - KLL = KLL.*polyval(LaguerrePoly(j),v_^2); - % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls - KdL1 = zeros(size(kx)); - for n_ = 0:Jmax-j - tmp_ = 0; - for s_ = 0:n_+j - tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2); - end - KdL1 = KdL1 + Kn(n_,kx).*tmp_; - end - % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls - KdL2 = zeros(size(kx)); - for n_ = 0:Jmax - tmp_ = 0; - for s_ = 0:min(Jmax,n_+j) - tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2); - end - KdL2 = KdL2 + Kn(n_,kx).*tmp_; - end - % errors - err_ = abs(J0Lj-KLL); - idx_ = 1; - while(err_(idx_)< eps && idx_ < numel(err_)) - idx_ = idx_ + 1; - end - klimLL(ij_,iv_) = kx(idx_); - err_ = abs(J0Lj-KdL1); - idx_ = 1; - while(err_(idx_)< eps && idx_ < numel(err_)) - idx_ = idx_ + 1; - end - klimdL1(ij_,iv_) = kx(idx_); - iv_ = iv_ + 1; - end -end -% plot -% plot(JmaxA,klimLL,'o-'); hold on -% plot(J_A,klimdL1,'o-'); hold on -end -surf(XX,YY,klimdL1) -xlabel('$J_{max}$'); ylabel('$v_{\perp}$'); -title(['$k$ s.t. $\epsilon=1\%$, $j=',num2str(j),'$']) -ylim([0,1]); grid on; -end if 0 %% Test Lj Ln = sum_s dnjs Ls -j = 10; -n = 10; -smax = n+j; -vperp = linspace(0,5,20); -LjLn = zeros(size(vperp)); -dnjsLs = zeros(size(vperp)); -dnjsLs_bin = zeros(size(vperp)); -dnjsLs_stir = zeros(size(vperp)); - -i=1; -for x_ = vperp - LjLn(i) = polyval(LaguerrePoly(j),x_)*polyval(LaguerrePoly(n),x_); - dnjsLs(i) = 0; - dnjsLs_bin(i) = 0; - dnjsLs_stir(i) = 0; +js = 4; +n = 4; +GYAC = 1; +smax = n+js; +wperp = linspace(0,5,50); +LjLn = zeros(size(wperp)); +dnjsLs = zeros(size(wperp)); +dnjsLs_bin = zeros(size(wperp)); +dnjsLs_stir = zeros(size(wperp)); +if GYAC + dnjsLs_GYAC = zeros(size(wperp)); + % Specify the filename + filename = '/home/ahoffman/gyacomo/results/dbg/output.txt'; + % Use the load function to read data from the text file + data = load(filename); + % Extract columns + indices = data(:, 1:3); + values = data(:, 4); + is_ = unique(indices(:,1)); + N_ = numel(is_); + dnjs_GYAC = zeros(N_,N_,N_); + for id2 = 1:numel(values) + in_ = indices(id2,1); + ij_ = indices(id2,2); + is_ = indices(id2,3); + dnjs_GYAC(in_,ij_,is_) = values(id2); + end +end + +id2=1; +for x_ = wperp + LjLn(id2) = polyval(LaguerrePoly(js),x_)*polyval(LaguerrePoly(n),x_); + dnjsLs(id2) = 0; + dnjsLs_bin(id2) = 0; + dnjsLs_stir(id2) = 0; for s_ = 0:smax - dnjsLs(i) = dnjsLs(i) + dnjs(n,j,s_)*polyval(LaguerrePoly(s_),x_); - dnjsLs_bin(i) = dnjsLs_bin(i) + dnjs_fact(n,j,s_)*polyval(LaguerrePoly(s_),x_); - dnjsLs_stir(i) = dnjsLs_stir(i) + dnjs_stir(n,j,s_)*polyval(LaguerrePoly(s_),x_); + dnjsLs(id2) = dnjsLs(id2) + dnjs(n,js,s_)*polyval(LaguerrePoly(s_),x_); + dnjsLs_bin(id2) = dnjsLs_bin(id2) + dnjs_fact(n,js,s_)*polyval(LaguerrePoly(s_),x_); + dnjsLs_stir(id2) = dnjsLs_stir(id2) + dnjs_stir(n,js,s_)*polyval(LaguerrePoly(s_),x_); + if GYAC + dnjsLs_GYAC(id2) = dnjsLs_GYAC(id2) + dnjs_GYAC(n+1,js+1,s_+1)*polyval(LaguerrePoly(s_),x_); + end end - i = i+1; + id2 = id2+1; end figure -plot(vperp,LjLn); hold on; -plot(vperp,dnjsLs,'d') -plot(vperp,dnjsLs_bin,'o') -plot(vperp,dnjsLs_stir/dnjsLs_stir(1),'x') +plot(wperp,LjLn,'DisplayName',['$L_',num2str(js),' L_',num2str(n),'$']); hold on; +plot(wperp,dnjsLs,'d','DisplayName',... + ['$\sum_{s=0}^{',num2str(smax),'} d^{coeff}_{',num2str(n),',',num2str(js),',s}L_s$']) +plot(wperp,dnjsLs_bin,'o','DisplayName','$\sum_s d^{fact}_{njs}L_s$') +plot(wperp,dnjsLs_stir/dnjsLs_stir(1),'x','DisplayName','$\sum_s d^{stir}_{njs}L_s$') +plot(wperp,dnjsLs_GYAC,'x','DisplayName','$GYAC$') +xlabel('$w_\perp$'); % We see that starting arround j = 18, n = 0, the stirling formula is the only % approximation that does not diverge from the target function. However it % performs badly for non 0 n... end - +if 0 +%% + wperp = linspace(0,3,128); + kp = linspace(0,3,128); + Jmax = 8; + js = Jmax/2;[0 2:Jmax]; + err_KLL = 1:numel(js); + err_KdL1 = 1:numel(js); + err_KdL2 = 1:numel(js); + err_T1 = 1:numel(js); + J0Lj = zeros(numel(kp),numel(wperp)); + KdL1 = zeros(numel(kp),numel(wperp)); + KdL2 = zeros(numel(kp),numel(wperp)); + T1 = zeros(numel(kp),numel(wperp)); + Kn = @(x__,y__) kernel(x__,y__); + % dnjs_ = @(n,j,s) dnjs(n,j,s); + dnjs_ = @(n,j,s) dnjs_GYAC(n+1,j+1,s+1); + for id1 = 1:numel(js) + for id2 = 1:numel(wperp) + j_ = js(id1); + v_ = wperp(id2); + % J0 x Lj + J0Lj_ = besselj(0,kp*v_)*polyval(LaguerrePoly(j_),v_^2); + % Taylor exp of J0 + T1_ = (1 - kp.^2*v_^2/4)*polyval(LaguerrePoly(j_),v_^2); + % Sum_n^Nmax Kn x Ln x Lj + KLL = zeros(size(kp)); + for n_ = 0:Jmax + KLL = KLL + Kn(n_,kp).*polyval(LaguerrePoly(n_),v_^2); + end + KLL = KLL.*polyval(LaguerrePoly(j_),v_^2); + % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls + KdL1_ = zeros(size(kp)); + for n_ = 0:Jmax-j_ + tmp_ = 0; + for s_ = 0:n_+j_ + tmp_ = tmp_ + dnjs_(n_,j_,s_)*polyval(LaguerrePoly(s_),v_^2); + end + KdL1_ = KdL1_ + Kn(n_,kp).*tmp_; + end + % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls + KdL2_ = zeros(size(kp)); + for n_ = 0:Jmax + tmp_ = 0; + for s_ = 0:min(Jmax,n_+j_) + tmp_ = tmp_ + dnjs_(n_,j_,s_)*polyval(LaguerrePoly(s_),v_^2); + end + KdL2_ = KdL2_ + Kn(n_,kp).*tmp_; + end + f_ = @(x) max(x); + kp2 = kp.^2+0.001; intJ0Lj = f_(abs(J0Lj_)); + err_KLL_(id2) = f_(abs(KLL -J0Lj_)./kp2)./f_(intJ0Lj./kp2); + err_KdL1_(id2) = f_(abs(KdL1_-J0Lj_)./kp2)./f_(intJ0Lj./kp2); + err_KdL2_(id2) = f_(abs(KdL2_-J0Lj_)./kp2)./f_(intJ0Lj./kp2); + err_T1_(id2) = f_(abs(T1_ -J0Lj_)./kp2)./f_(intJ0Lj./kp2); + J0Lj(:,id2) = J0Lj_; + KdL1(:,id2) = KdL1_; + KdL2(:,id2) = KdL2_; + T1 (:,id2) = T1_; + end + err_KLL(id1) = sum(err_KLL_)./numel(wperp); + err_KdL1(id1) = sum(err_KdL1_)./numel(wperp); + err_KdL2(id1) = sum(err_KdL2_)./numel(wperp); + err_T1(id1) = sum(err_T1_)./numel(wperp); + end + clr_ = lines(10); + figure + plot(js,err_KLL,'-',... + 'DisplayName','$\sum_{n=0}^{J}\mathcal K_n L_n L_j$'); hold on + plot(js,err_KdL1,':x',... + 'DisplayName','$\sum_{n=0}^{J-j}\mathcal K_n \sum_{s=0}^{n+j} d_{njs} L_s$',... + 'Color',clr_(5,:)); hold on + plot(js,err_KdL2,':o',... + 'DisplayName','$\sum_{n=0}^{J}\mathcal K_n \sum_{s=0}^{\min(J,n+j)} d_{njs} L_s$',... + 'Color',clr_(2,:)); hold on + plot(js,err_T1,':s',... + 'DisplayName','$(1 - l_\perp) L_j$',... + 'Color',clr_(4,:)); hold on + set(gca,'YScale','log'); + xlabel('$j$'); ylabel('$\varepsilon_J$') + if 1 + %% + figure + nc = 25 + subplot(131) + contourf(kp,wperp,J0Lj',nc); xlabel('$l_{\perp a}$'); ylabel('$w_{\perp a}$') + clim_ = clim; + subplot(132) + contourf(kp,wperp,KdL1', nc); xlabel('$l_{\perp a}$'); ylabel('$w_{\perp a}$') + clim(clim_); + % subplot(132) + subplot(133) + contourf(kp,wperp,KdL2', nc); xlabel('$l_{\perp a}$'); ylabel('$w_{\perp a}$') + % contourf(kp,wperp,T1 ,20); xlabel('$k_\perp$'); ylabel('$w_\perp$') + clim(clim_); + end +end diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m index 1267999b25d0121511e955858f9b6b53c6509e4f..9a88d3b5dae77018ff8478197e9bd7e90b43ef9b 100644 --- a/matlab/plot/plot_metric.m +++ b/matlab/plot/plot_metric.m @@ -96,9 +96,14 @@ if NPLOT > 0 end end %outputs -arrays = squeeze(geo_arrays(:,:)); +% arrays = squeeze(geo_arrays(:,:)); R = [geo_arrays(:,12);geo_arrays(1,12)]; Z = [geo_arrays(:,13);geo_arrays(1,13)]; + +for i_ = 1:numel(names) + namae = names{i_}; + arrays.(namae) = geo_arrays(:,i_); +end try if ~options.SHOWFIG close diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m index 8fdcbc3d2cb4250cb560cd2de67a429ea9032e70..a33d3612a46247dd21d15ac9e292f418a18fbbe8 100644 --- a/matlab/plot/plot_radial_transport_and_spacetime.m +++ b/matlab/plot/plot_radial_transport_and_spacetime.m @@ -56,7 +56,8 @@ mvm = @(x) movmean(x,OPTIONS.NMVA); grid on; set(gca,'xticklabel',[]); title({DATA.param_title,... ['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]}); -%% radial shear radial profile + axis tight + %% radial shear radial profile % computation Ns3D = numel(DATA.Ts3D); [KX, KY] = meshgrid(DATA.grids.kx, DATA.grids.ky); @@ -68,18 +69,33 @@ mvm = @(x) movmean(x,OPTIONS.NMVA); OPTIONS.NAME = OPTIONS.ST_FIELD; OPTIONS.PLAN = 'xy'; OPTIONS.COMP = 'avg'; + OPTIONS.TAVG = 0; OPTIONS.TIME = DATA.Ts3D; OPTIONS.POLARPLOT = 0; - toplot = process_field(DATA,OPTIONS); - f2plot = toplot.FIELD; + try + if DATA.fort_00.DIAGNOSTICS.dtsave_2d > 0 + opt_.PLAY = 0; + [MOVIE] = cinecita_2D(DATA.dir,OPTIONS.ST_FIELD,opt_); + f2plot = MOVIE.F; + TIME = MOVIE.T; + else + toplot = process_field(DATA,OPTIONS); + f2plot = toplot.FIELD; + TIME = DATA.Ts3D(toplot.FRAMES); + end + catch + toplot = process_field(DATA,OPTIONS); + f2plot = toplot.FIELD; + TIME = DATA.Ts3D(toplot.FRAMES); + end dframe = ite3D - its3D; - clim = max(max(max(abs(plt(f2plot(:,:,:)))))); + clim_ = max(max(max(abs(plt(f2plot(:,:,:)))))); FIGURE.ax3 = subplot(2,1,2,'parent',FIGURE.fig); - [TY,TX] = meshgrid(DATA.grids.x,DATA.Ts3D(toplot.FRAMES)); + [TY,TX] = meshgrid(DATA.grids.x,TIME); pclr = pcolor(TX,TY,squeeze(plt(f2plot))'); set(pclr, 'edgecolor','none'); legend(['$\langle ',OPTIONS.ST_FIELD,' \rangle_{y,z}$']) %colorbar; - caxis(clim*[-1 1]); + clim(clim_*[-1 1]); cmap = bluewhitered(256); colormap(cmap) xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m index ab9d3aea55bd247e410291c6560aa38f823805d4..f193f118b32725c1187cee06d4fe97f8b1bf215e 100644 --- a/matlab/plot/show_moments_spectrum.m +++ b/matlab/plot/show_moments_spectrum.m @@ -6,16 +6,8 @@ Ja = DATA.grids.Jarray; Time_ = DATA.Ts3D; FIGURE.fig = figure; FIGURE.FIGNAME = ['mom_spectrum_',DATA.params_string]; % set(gcf, 'Position', [100 50 1000 400]) -if OPTIONS.ST == 0 - OPTIONS.LOGSCALE = 0; -end -if OPTIONS.LOGSCALE - logname = 'log'; - compress = @(x,ia) log(sum(abs((x(ia,:,:,:,:))),4)); -else - logname = ''; - compress = @(x,ia) sum(abs((x(ia,:,:,:,:))),4); -end +logname = ''; +compress = @(x,ia) sum(abs((x(ia,:,:,:,:))),4); for ia = 1:DATA.inputs.Na Napjz = compress(DATA.Napjz,ia); Napjz = reshape(Napjz,DATA.grids.Np,DATA.grids.Nj,numel(DATA.Ts3D)); @@ -71,8 +63,10 @@ for ia = 1:DATA.inputs.Na % plots % scaling if OPTIONS.TAVG_2D - nt_half = ceil(numel(Time_)/2); - Napjz_tavg = mean(Napjz(:,:,nt_half:end),3); + t0 = OPTIONS.TWINDOW(1); t1 = OPTIONS.TWINDOW(2); + [~,it0] = min(abs(DATA.Ts3D-t0)); + [~,it1] = min(abs(DATA.Ts3D-t1)); + Napjz_tavg = mean(Napjz(:,:,it0:it1),3); x_ = DATA.grids.Parray; y_ = DATA.grids.Jarray; if OPTIONS.TAVG_2D_CTR @@ -84,7 +78,7 @@ for ia = 1:DATA.inputs.Na xlabel('$p$'); ylabel('$j$'); end clb = colorbar; colormap(jet); - clb.Label.String = '$\log\langle E(p,j) \rangle_t$'; + clb.Label.String = '$\langle E(p,j) \rangle_t$'; clb.TickLabelInterpreter = 'latex'; clb.Label.Interpreter = 'latex'; clb.Label.FontSize= 18; @@ -95,7 +89,7 @@ for ia = 1:DATA.inputs.Na plt = @(x,i) reshape(x(i,:),numel(p2j),numel(Time_)); end Nlabelmax = 15; - nskip = floor(numel(p2j)/Nlabelmax); + nskip = max(1,floor(numel(p2j)/Nlabelmax)); if OPTIONS.ST imagesc(Time_,p2j,plt(Na_ST,1:numel(p2j))); set(gca,'YDir','normal') @@ -105,10 +99,13 @@ for ia = 1:DATA.inputs.Na yticklabels(ticks_labels(1:nskip:end)) end if OPTIONS.FILTER - caxis([0 0.2]); - title('Relative fluctuation from the average'); + clim([0 0.2]); + title('Relative fluctuation from the average'); end colorbar + if OPTIONS.LOGSCALE + set(gca,'ColorScale','log'); + end else colors_ = jet(numel(p2j)); for i = 1:numel(p2j) @@ -120,6 +117,7 @@ for ia = 1:DATA.inputs.Na end end top_title(DATA.paramshort) + title(plotname) end diff --git a/matlab/plot/tight_layout_generator.m b/matlab/plot/tight_layout_generator.m new file mode 100644 index 0000000000000000000000000000000000000000..f85c09ea5da8cb84dd92c719614570d0a22ecdf7 --- /dev/null +++ b/matlab/plot/tight_layout_generator.m @@ -0,0 +1,11 @@ +function [] = tight_layout_generator(m,n) +% Generate a layout with tight subplots +% m = 3; % number lines +% n = 3; % number of columns +figure +t = tiledlayout(m,n,'TileSpacing','Compact','Padding','Compact'); +for i = 1:(m*n) +nexttile +end + +end diff --git a/matlab/process_field.m b/matlab/process_field.m index 9390b1b2cd33e40186f458a292373df32dad64bd..9ead884c631372d4790cf003240b3e5aebbad2df 100644 --- a/matlab/process_field.m +++ b/matlab/process_field.m @@ -8,45 +8,44 @@ for i = 1:numel(OPTIONS.TIME) [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D)); end FRAMES = unique(FRAMES); +TIME = DATA.Ts3D(FRAMES); %% Setup the plot geometry [KX, KY] = meshgrid(DATA.grids.kx, DATA.grids.ky); directions = {'y','x','z'}; Nx = DATA.grids.Nx; Ny = DATA.grids.Ny; Nz = DATA.grids.Nz; Nt = numel(FRAMES); -POLARPLOT = OPTIONS.POLARPLOT; LTXNAME = OPTIONS.NAME; +SKIP_COMP = 0; % Turned on only for kin. distr. func. plot; +Z = 0; switch OPTIONS.PLAN case 'xy' - XNAME = '$x$'; YNAME = '$y$'; + XNAME = '$x/\rho_s$'; YNAME = '$y/\rho_s$'; [X,Y] = meshgrid(DATA.grids.x,DATA.grids.y); REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1; case 'xz' - XNAME = '$x$'; YNAME = '$z$'; + XNAME = '$x/\rho_s$'; YNAME = '$z/L_\parallel$'; [Y,X] = meshgrid(DATA.grids.z,DATA.grids.x); - REALP = 1; COMPDIM = 1; SCALE = 0; + REALP = 1; COMPDIM = 1; SCALE = 0; POLARPLOT = 0; case 'yz' - XNAME = '$y$'; YNAME = '$z$'; + XNAME = '$y/\rho_s$'; YNAME = '$z/L_\parallel$'; [Y,X] = meshgrid(DATA.grids.z,DATA.grids.y); - REALP = 1; COMPDIM = 2; SCALE = 0; + REALP = 1; COMPDIM = 2; POLARPLOT = 0; SCALE = 0; case 'kxky' - XNAME = '$k_x$'; YNAME = '$k_y$'; + XNAME = '$k_x\rho_s$'; YNAME = '$k_y\rho_s$'; [X,Y] = meshgrid(DATA.grids.kx,DATA.grids.ky); REALP = 0; COMPDIM = 3; POLARPLOT = 0; SCALE = 1; case 'kxz' - XNAME = '$k_x$'; YNAME = '$z$'; + XNAME = '$k_x\rho_s$'; YNAME = '$z/L_\parallel$'; [Y,X] = meshgrid(DATA.grids.z,DATA.grids.kx); REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0; case 'kyz' - XNAME = '$k_y$'; YNAME = '$z$'; + XNAME = '$k_y\rho_s$'; YNAME = '$z/L_\parallel$'; [Y,X] = meshgrid(DATA.grids.z,DATA.grids.ky); REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0; - case 'sx' - XNAME = '$v_\parallel$'; YNAME = '$\mu$'; - [Y,X] = meshgrid(OPTIONS.XPERP,OPTIONS.SPAR); - REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 0; - for i = 1:numel(OPTIONS.TIME) - [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts5D)); - end - FRAMES = unique(FRAMES); Nt = numel(FRAMES); + case 'xyz' + XNAME = '$x/\rho_s$'; YNAME = '$y/\rho_s$'; + [X,Y] = meshgrid(DATA.grids.x,DATA.grids.y); + REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1; + SKIP_COMP = 1; OPTIONS.COMP = 3; end Xmax_ = max(max(abs(X))); Ymax_ = max(max(abs(Y))); Lmin_ = min([Xmax_,Ymax_]); @@ -73,29 +72,6 @@ else FIELD = zeros(sz(1),sz(2),Nt); end %% Process the field to plot -% short term writing -- -% b_e = DATA.sigma_e*sqrt(2*DATA.tau_e)*sqrt(KX.^2+KY.^2); -% adiab_e = kernel(0,b_e); -% pol_e = kernel(0,b_e).^2; -% for n = 1:DATA.Jmaxe -% adiab_e = adiab_e + kernel(n,b_e).^2; -% pol_e = pol_e + kernel(n,b_e).^2; -% end -% adiab_e = DATA.q_e/DATA.tau_e .* adiab_e; -% pol_e = DATA.q_e^2/DATA.tau_e * (1 - pol_e); -% -% b_i = DATA.sigma_i*sqrt(2*DATA.tau_i)*sqrt(KX.^2+KY.^2); -% adiab_i = kernel(0,b_i); -% pol_i = kernel(0,b_i).^2; -% for n = 1:DATA.Jmaxi -% adiab_i = adiab_i + kernel(n,b_i).^2; -% pol_i = pol_i + kernel(n,b_i).^2; -% end -% pol_i = DATA.q_i^2/DATA.tau_i * (1 - pol_i); -% adiab_i = DATA.q_i/DATA.tau_i .* adiab_i; -% poisson_op = (pol_i + pol_e); -adiab_e =0; adiab_i =0; poisson_op=1; -SKIP_COMP = 0; % Turned on only for kin. distr. func. plot % -- switch OPTIONS.COMP case 'avg' @@ -112,18 +88,18 @@ switch OPTIONS.COMP i = OPTIONS.COMP; compr = @(x) x(i,:,:); if REALP - COMPNAME = sprintf(['y=','%2.1f'],DATA.grids.x(i)); + COMPNAME = sprintf(['y=','%2.1f'],DATA.grids.y(i)); else - COMPNAME = sprintf(['k_y=','%2.1f'],DATA.grids.kx(i)); + COMPNAME = sprintf(['k_y=','%2.1f'],DATA.grids.ky(i)); end FIELDNAME = [LTXNAME,'(',COMPNAME,')']; case 2 i = OPTIONS.COMP; compr = @(x) x(:,i,:); if REALP - COMPNAME = sprintf(['x=','%2.1f'],DATA.grids.y(i)); + COMPNAME = sprintf(['x=','%2.1f'],DATA.grids.x(i)); else - COMPNAME = sprintf(['k_x=','%2.1f'],DATA.grids.ky(i)); + COMPNAME = sprintf(['k_x=','%2.1f'],DATA.grids.kx(i)); end FIELDNAME = [LTXNAME,'(',COMPNAME,')']; case 3 @@ -188,6 +164,10 @@ switch OPTIONS.NAME NAME = 'sx'; FLD_ = DATA.PHI(:,:,:,FRAMES); OPE_ = KY.^2; + case 'w_{Ez}' % x-comp of ExB shear + NAME = 'wEz'; + FLD_ = DATA.PHI(:,:,:,FRAMES); + OPE_ = (KX.^2+KY.^2).*(abs(KX)<=2).*(abs(KY)<=2); case '\omega_z' % ES pot vorticity NAME = 'vorticity'; FLD_ = DATA.PHI(:,:,:,FRAMES); @@ -206,11 +186,11 @@ switch OPTIONS.NAME OPE_ = 1; case 'n_i' % ion density NAME = 'ni'; - FLD_ = DATA.DENS_I(:,:,:,FRAMES) - adiab_i.* DATA.PHI(:,:,:,FRAMES); + FLD_ = DATA.DENS_I(:,:,:,FRAMES); OPE_ = 1; case 'n_e' % electron density NAME = 'ne'; - FLD_ = DATA.DENS_E(:,:,:,FRAMES) - adiab_e.* DATA.PHI(:,:,:,FRAMES); + FLD_ = DATA.DENS_E(:,:,:,FRAMES); OPE_ = 1; case 'k^2n_e' % electron vorticity NAME = 'k2ne'; @@ -219,12 +199,24 @@ switch OPTIONS.NAME case 'n_i-n_e' % polarisation NAME = 'pol'; OPE_ = 1; - FLD_ = ((DATA.DENS_I(:,:,:,FRAMES)- adiab_i.* DATA.PHI(:,:,:,FRAMES))... - -(DATA.DENS_E(:,:,:,FRAMES)- adiab_e.* DATA.PHI(:,:,:,FRAMES))); + FLD_ = DATA.DENS_I(:,:,:,FRAMES)... + -DATA.DENS_E(:,:,:,FRAMES); + case 'u_i' % ion density + NAME = 'ui'; + FLD_ = DATA.UPAR_I(:,:,:,FRAMES); + OPE_ = 1; + case 'u_e' % electron density + NAME = 'ue'; + FLD_ = DATA.UPAR_E(:,:,:,FRAMES); + OPE_ = 1; case 'T_i' % ion temperature NAME = 'Ti'; FLD_ = DATA.TEMP_I(:,:,:,FRAMES); OPE_ = 1; + case 'T_e' % ion temperature + NAME = 'Te'; + FLD_ = DATA.TEMP_E(:,:,:,FRAMES); + OPE_ = 1; case 'G_x' % ion particle flux NAME = 'Gx'; FLD_ = 0.*DATA.PHI(:,:,:,FRAMES); @@ -238,24 +230,55 @@ switch OPTIONS.NAME % tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.grids.Ny,DATA.grids.Nx))),2)); tmp(:,:,iz) = abs((squeeze(fft2(gx_,DATA.grids.Ny,DATA.grids.Nx)))); end - FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:))); - end - case 'Q_x' % ion heat flux - NAME = 'Qx'; - % FLD_ = 0.*DATA.PHI(:,:,:,FRAMES); + FLD_(:,:,:,it)= squeeze((tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:))); + end + case 'n_i T_i' % ion heat flux + NAME = 'niTi'; + FLD_ = 0.*DATA.PHI(:,:,:,FRAMES); OPE_ = 1; for it = 1:numel(FRAMES) - tmp = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz); + qx_c = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz); for iz = 1:DATA.grids.Nz - vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it)))))); - ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it))))); - Ti_ = real((ifourier_GENE( DATA.TEMP_I(:,:,iz,FRAMES(it))))); - qx_ = vx_.*ni_.*Ti_; -% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.grids.Ny,DATA.grids.Nx))),2)); - tmp(:,:,iz) = abs((squeeze(fft2(qx_,DATA.grids.Ny,DATA.grids.Nx)))); + ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it))))); + Ti_ = real((ifourier_GENE( DATA.TEMP_I(:,:,iz,FRAMES(it))))); + cx_r = ni_.*Ti_; + cx_ = fft(cx_r,[],1); + cx_c(:,:,iz) = fft(cx_ ,[],2); end - FLD_(:,:,:,it)= squeeze(tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:)); - end + FLD_(:,:,:,it)= squeeze(cx_c(1:DATA.grids.Nky,1:DATA.grids.Nkx,:))./numel(qx_c); + end + case 'Q_{xi}' % ion heat flux + NAME = 'Qxi'; + FLD_ = 0.*DATA.PHI(:,:,:,FRAMES); + OPE_ = 1; + for it = 1:numel(FRAMES) + qx_c = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz); + for iz = 1:DATA.grids.Nz + vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it)))))); + ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it))))); + Ti_ = real((ifourier_GENE( DATA.TEMP_I(:,:,iz,FRAMES(it))))); + qx_r = vx_.*ni_.*Ti_; + qx_ = fft(qx_r,[],1); + qx_c(:,:,iz) = fft(qx_ ,[],2); + end + FLD_(:,:,:,it)= squeeze(qx_c(1:DATA.grids.Nky,1:DATA.grids.Nkx,:))./numel(qx_c); + end + case 'Q_{xe}' % electron heat flux + NAME = 'Qxe'; + FLD_ = 0.*DATA.PHI(:,:,:,FRAMES); + OPE_ = 1; + for it = 1:numel(FRAMES) + qx_c = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz); + for iz = 1:DATA.grids.Nz + vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it)))))); + ne_ = real((ifourier_GENE( DATA.DENS_E(:,:,iz,FRAMES(it))))); + Te_ = real((ifourier_GENE( DATA.TEMP_E(:,:,iz,FRAMES(it))))); + qx_r = vx_.*ne_.*Te_; + qx_ = fft(qx_r,[],1); + qx_c(:,:,iz) = fft(qx_ ,[],2); + end + FLD_(:,:,:,it)= squeeze(qx_c(1:DATA.grids.Nky,1:DATA.grids.Nkx,:))./numel(qx_c); + end case 'f_i' SKIP_COMP = 1; shift_x = @(x) x; shift_y = @(y) y; @@ -284,29 +307,34 @@ switch OPTIONS.NAME end % Process the field according to the 2D plane and the space (real/cpx) if ~SKIP_COMP -if COMPDIM == 3 - for it = 1:numel(FRAMES) - tmp = squeeze(compr(OPE_.*FLD_(:,:,:,it))); - FIELD(:,:,it) = squeeze(process(tmp)); - end -else - if REALP - tmp = zeros(Ny,Nx,Nz); + if COMPDIM == 3 + for it = 1:numel(FRAMES) + tmp = squeeze(compr(OPE_.*FLD_(:,:,:,it))); + FIELD(:,:,it) = squeeze(process(tmp)); + end else - tmp = zeros(DATA.grids.Nky,DATA.grids.Nkx,Nz); - end - for it = 1:numel(FRAMES) - for iz = 1:numel(DATA.grids.z) - tmp(:,:,iz) = squeeze(process(OPE_.*FLD_(:,:,iz,it))); + if REALP + tmp = zeros(Ny,Nx,Nz); + else + tmp = zeros(DATA.grids.Nky,DATA.grids.Nkx,Nz); end - FIELD(:,:,it) = squeeze(compr(tmp)); - end -end + for it = 1:numel(FRAMES) + for iz = 1:numel(DATA.grids.z) + tmp(:,:,iz) = squeeze(process(OPE_.*FLD_(:,:,iz,it))); + end + FIELD(:,:,it) = squeeze(compr(tmp)); + end + end +else + clear FIELD; + FIELD = squeeze(process(FLD_)); end TOPLOT.FIELD = FIELD; +TOPLOT.TIME = TIME; TOPLOT.FRAMES = FRAMES; TOPLOT.X = shift_x(X); TOPLOT.Y = shift_y(Y); +TOPLOT.Z = Z; TOPLOT.FIELDNAME = FIELDNAME; TOPLOT.XNAME = XNAME; TOPLOT.YNAME = YNAME; @@ -315,5 +343,6 @@ TOPLOT.DIMENSIONS= DIMENSIONS; TOPLOT.ASPECT = ASPECT; TOPLOT.FRAMES = FRAMES; TOPLOT.INTERP = INTERP; + end diff --git a/matlab/read_flux_out_XX.m b/matlab/read_flux_out_XX.m index 2c7737ea4c28c895ca065fcae4df59af0b02c592..22f4ce581ac07712ec5492602e0ff5f354c28395 100644 --- a/matlab/read_flux_out_XX.m +++ b/matlab/read_flux_out_XX.m @@ -1,4 +1,4 @@ -function [t_all, Pxi_all, Qxi_all, Pxe_all, Qxe_all] = read_flux_out_XX(folderPath,PLOT) +function [out] = read_flux_out_XX(folderPath,PLOT) % Check if the prompt_string is provided as an argument if nargin < 1 % If not provided, prompt the user for input @@ -90,4 +90,11 @@ function [t_all, Pxi_all, Qxi_all, Pxe_all, Qxe_all] = read_flux_out_XX(folderPa legend('show'); title('Radial heat flux') end + + + out.t = t_all; + out.Pxi = Pxi_all; + out.Qxi = Qxi_all; + out.Pxe = Pxe_all; + out.Qxe = Qxe_all; end diff --git a/src/coeff_mod.F90 b/src/coeff_mod.F90 index 3c4c9125a19bfb7103e6050140f4850558bb715e..9739766ce2f727db1bac0fc445abd99484458646 100644 --- a/src/coeff_mod.F90 +++ b/src/coeff_mod.F90 @@ -1,7 +1,9 @@ MODULE coeff -! this module contains routines to compute normalization coefficients and velocity integrals -! +! this module contains routines to compute the Laguerre-Laguerre products coefficients +! dnjs for the nonlinear term. ALL2L(n,j,s,0) = dnjs +! Lpjl(p,j,l) returns the l-th coefficient of the j-th order associated Laguerre_{p-1/2} poynomial + ! the canonical basis (p=0) writes L_j(x) = sum_{l=0}^j c_l x^l USE PREC_CONST use BASIC diff --git a/src/initial_mod.F90 b/src/initial_mod.F90 index a744d64ee3fd253ae3ab52d8dad7af984cb4e9da..94eaedb74026ce00db61b6bb5869b109be7a44b0 100644 --- a/src/initial_mod.F90 +++ b/src/initial_mod.F90 @@ -289,8 +289,8 @@ CONTAINS USE fields, ONLY: moments USE prec_const, ONLY: xp USE model, ONLY: LINEARITY - USE grid, ONLY: total_nkx, local_nkx_offset, local_nky, local_nky_offset,& - kxarray_full, kyarray_full, kx_max, kx_min, ky_max + USE grid, ONLY: total_nkx, local_nkx_offset, local_nky, local_nky_offset, iky0,& + kxarray_full, kyarray_full, kx_max, kx_min, ky_max, deltakx, deltaky IMPLICIT NONE INTEGER :: ikx,iky, im, I_, J_ REAL(xp):: maxkx, minkx, maxky, kx_, ky_, A_ @@ -314,19 +314,16 @@ CONTAINS IF(I_ .LT. 0) THEN I_ = I_ + total_nkx ENDIF - I_ = I_ + 1 ! array indices start at 1 - J_ = J_ + 1 ! Check the validity of the modes - kx_ = kxarray_full(I_) - ky_ = kyarray_full(J_) + kx_ = I_*deltakx + ky_ = J_*deltaky IF( ((kx_ .GE. minkx) .AND. (kx_ .LE. maxkx)) .AND. & ((ky_ .GE. 0._xp) .AND. (ky_ .LE. maxky)) ) THEN ! Init the mode DO ikx=1,total_nkx DO iky=1,local_nky - IF ( (ikx+local_nkx_offset .EQ. I_) .AND. & - (iky+local_nky_offset .EQ. J_) ) THEN - ! WRITE(*,'(A10,F4.2,A,F4.2,A,F4.2)') '-init (kx=',kx_,',ky=',ky_,') at Amp=',A_ + IF ( (kxarray_full(ikx+local_nkx_offset) .EQ. kx_) .AND. & + (kyarray_full(iky+local_nky_offset) .EQ. ky_) ) THEN WRITE(*,'(A,F5.3,A,F5.3,A,G9.2)') '-init (kx=',kx_,',ky=',ky_,') with Amp= ',A_ moments(:,:,:,iky,ikx,:,:) = A_ ENDIF diff --git a/src/model_mod.F90 b/src/model_mod.F90 index 81ebdbe7494f15083d894913421ed1ef707d8554..7edf6e87ebaebf62ab330bd3003edbbe494bc1d6 100644 --- a/src/model_mod.F90 +++ b/src/model_mod.F90 @@ -4,7 +4,6 @@ MODULE model IMPLICIT NONE PRIVATE ! INPUTS - INTEGER, PUBLIC, PROTECTED :: KERN = 0 ! Kernel model CHARACTER(len=16), & PUBLIC, PROTECTED ::LINEARITY= 'linear' ! To turn on non linear bracket term REAL(xp), PUBLIC, PROTECTED :: mu_x = 0._xp ! spatial x-Hyperdiffusivity coefficient (for num. stability) @@ -39,6 +38,11 @@ MODULE model LOGICAL, PUBLIC, PROTECTED :: FORCE_SYMMETRY = .false. ! Add or remove the ExB nonlinear correction (Mcmillan 2019) LOGICAL, PUBLIC, PROTECTED :: ExB_NL_CORRECTION = .true. + CHARACTER(len=16), & + PUBLIC, PROTECTED :: KN_MODEL = 'std' ! Kernel model + INTEGER, PUBLIC, PROTECTED :: ORDER = 1 ! order for Taylor expansion + INTEGER, PUBLIC, PROTECTED :: ORDER_NUM = 2 ! numerator order for Pade approx + INTEGER, PUBLIC, PROTECTED :: ORDER_DEN = 4 ! denominator order for Pade approx ! Module's routines PUBLIC :: model_readinputs, model_outputinputs @@ -52,11 +56,11 @@ CONTAINS USE prec_const, ONLY: xp IMPLICIT NONE - NAMELIST /MODEL/ KERN, LINEARITY, RM_LD_T_EQ, FORCE_SYMMETRY, MHD_PD, & + NAMELIST /MODEL/ LINEARITY, RM_LD_T_EQ, FORCE_SYMMETRY, MHD_PD, & Na, ADIAB_E, ADIAB_I, q_i, tau_i, & mu_x, mu_y, N_HD, HDz_h, mu_z, mu_p, mu_j, HYP_V, & nu, k_gB, k_cB, lambdaD, beta, ExBrate, ExB_NL_CORRECTION,& - ikxZF, ZFrate, ZF_ONLY + ikxZF, ZFrate, ZF_ONLY, KN_MODEL, ORDER, ORDER_NUM, ORDER_DEN READ(lu_in,model) @@ -96,7 +100,6 @@ CONTAINS CHARACTER(len=256) :: str WRITE(str,'(a)') '/data/input/model' CALL creatd(fid, 0,(/0/),TRIM(str),'Model Input') - CALL attach(fid, TRIM(str), "KERN", KERN) CALL attach(fid, TRIM(str), "LINEARITY", LINEARITY) CALL attach(fid, TRIM(str),"RM_LD_T_EQ",RM_LD_T_EQ) CALL attach(fid, TRIM(str), "mu_x", mu_x) @@ -120,6 +123,10 @@ CONTAINS CALL attach(fid, TRIM(str), "ADIAB_E", ADIAB_E) CALL attach(fid, TRIM(str), "ADIAB_I", ADIAB_I) CALL attach(fid, TRIM(str), "tau_i", tau_i) + CALL attach(fid, TRIM(str),"kern model",KN_MODEL) + CALL attach(fid, TRIM(str), "order den",ORDER_DEN) + CALL attach(fid, TRIM(str), "order num",ORDER_NUM) + CALL attach(fid, TRIM(str), "order", ORDER) END SUBROUTINE model_outputinputs END MODULE model diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90 index c11fa8f8422c98b2726d56399645652cf496d016..7b8fb4308df6015832cb7c77caa521821caa413c 100644 --- a/src/nonlinear_mod.F90 +++ b/src/nonlinear_mod.F90 @@ -94,6 +94,7 @@ SUBROUTINE compute_nonlinear ENDIF ! Second convolution terms G_cmpx = 0._xp ! initialization of the sum + ! We set smax as min(n+j,Jmax) for safety reason but nmaxarray may prevent it as well (see anti-Laguerre-aliasing truncation) smax = MIN( jarray(ini)+jarray(iji), jmax ); s1:DO is = 1, smax+1 ! sum truncation on number of moments isi = is + ngj/2 diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90 index 62f242199ff654f81ec877bf498c3a6edf393283..357f092c1bdd1222e59aa38b6b6b13f4ca6ace2e 100644 --- a/src/numerics_mod.F90 +++ b/src/numerics_mod.F90 @@ -77,31 +77,80 @@ SUBROUTINE evaluate_kernels USE grid, ONLY : local_na, local_nj,ngj, local_nkx, local_nky, local_nz, ngz, jarray, kp2array,& nzgrid USE species, ONLY : sigma2_tau_o2 + USE model, ONLY : KN_MODEL, ORDER, ORDER_NUM, ORDER_DEN IMPLICIT NONE INTEGER :: j_int, ia, eo, ikx, iky, iz, ij REAL(xp) :: j_xp, y_, factj, sigma_i sigma_i = 1._xp ! trivial singe sigma_a = sqrt(m_a/m_i) -DO ia = 1,local_na - DO eo = 1,nzgrid - DO ikx = 1,local_nkx - DO iky = 1,local_nky - DO iz = 1,local_nz + ngz - DO ij = 1,local_nj + ngj - j_int = jarray(ij) - j_xp = REAL(j_int,xp) - y_ = sigma2_tau_o2(ia) * kp2array(iky,ikx,iz,eo) - IF(j_int .LT. 0) THEN !ghosts values - kernel(ia,ij,iky,ikx,iz,eo) = 0._xp - ELSE - factj = REAL(GAMMA(j_xp+1._xp),xp) - kernel(ia,ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/factj - ENDIF + SELECT CASE (KN_MODEL) + CASE('taylor') ! developped with Leonhard Driever + ! Kernels based on the ORDER_NUM order Taylor series of J0 + WRITE (*,*) 'Kernel approximation uses Taylor series with ', ORDER, ' powers of k' + DO ia = 1,local_na + DO eo = 1,nzgrid + DO ikx = 1,local_nkx + DO iky = 1,local_nky + DO iz = 1,local_nz + ngz + DO ij = 1,local_nj + ngj + y_ = sigma2_tau_o2(ia) * kp2array(iky,ikx,iz,eo) + j_int = jarray(ij) + IF (j_int > ORDER .OR. j_int < 0) THEN + kernel(ia,ij,ikx,iky,iz,eo) = 0._xp + ELSE + kernel(ia,ij,ikx,iky,iz,eo) = taylor_kernel_n(ORDER, j_int, y_) + ENDIF + ENDDO + ENDDO ENDDO ENDDO ENDDO ENDDO - ENDDO + CASE ('pade') + ! Kernels based on the ORDER_NUM / ORDER_DEN Pade approximation of the kernels + WRITE (*,*) 'Kernel approximation uses ', ORDER_NUM ,'/', ORDER_DEN, ' Pade approximation' + DO ia = 1,local_na + DO eo = 1,nzgrid + DO ikx = 1,local_nkx + DO iky = 1,local_nky + DO iz = 1,local_nz + ngz + DO ij = 1,local_nj + ngj + y_ = sigma2_tau_o2(ia) * kp2array(iky,ikx,iz,eo) + j_int = jarray(ij) + IF (j_int > ORDER_NUM .OR. j_int < 0) THEN + kernel(ia,ij,ikx,iky,iz,eo) = 0._xp + ELSE + kernel(ia,ij,ikx,iky,iz,eo) = pade_kernel_n(j_int, y_,ORDER_NUM,ORDER_DEN) + ENDIF + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + CASE DEFAULT + DO ia = 1,local_na + DO eo = 1,nzgrid + DO ikx = 1,local_nkx + DO iky = 1,local_nky + DO iz = 1,local_nz + ngz + DO ij = 1,local_nj + ngj + j_int = jarray(ij) + j_xp = REAL(j_int,xp) + y_ = sigma2_tau_o2(ia) * kp2array(iky,ikx,iz,eo) + IF(j_int .LT. 0) THEN !ghosts values + kernel(ia,ij,iky,ikx,iz,eo) = 0._xp + ELSE + factj = REAL(GAMMA(j_xp+1._xp),xp) + kernel(ia,ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/factj + ENDIF + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO + END SELECT ! !! Correction term for the evaluation of the heat flux ! HF_phi_correction_operator(:,:,:) = & @@ -117,7 +166,6 @@ DO ia = 1,local_na ! - (j_xp+1.0_xp) * Kernel(ia,ij+1,:,:,:,1) & ! - j_xp * Kernel(ia,ij-1,:,:,:,1)) ! ENDDO -ENDDO END SUBROUTINE evaluate_kernels !******************************************************************************! @@ -324,4 +372,123 @@ SUBROUTINE compute_lin_coeff ENDDO END SUBROUTINE compute_lin_coeff + +!******************************************************************************! +!!!!!!! Auxilliary kernel functions/subroutines (developped with Leonhard Driever) +!******************************************************************************! +REAL(xp) FUNCTION taylor_kernel_n(order, n, y) + IMPLICIT NONE + INTEGER, INTENT(IN) :: order + INTEGER, INTENT(IN) :: n + REAL(xp), INTENT(IN) :: y + REAL(xp) :: sum_variable + INTEGER :: m + REAL(xp) :: m_dp, n_dp + + n_dp = REAL(n, xp) + sum_variable = 0._xp + + DO m = n, order + m_dp = REAL(m, xp) + sum_variable = sum_variable + (-1._xp)**(m - n) * y**m / (GAMMA(n_dp + 1._xp) * GAMMA(m_dp - n_dp + 1._xp)) ! Denominator of m C n + END DO + taylor_kernel_n = sum_variable +END FUNCTION taylor_kernel_n + + +REAL(xp) FUNCTION pade_kernel_n(n, y, N_NUM, N_DEN) + IMPLICIT NONE + INTEGER, INTENT(IN) :: n, N_NUM, N_DEN + REAL(xp), INTENT(IN) :: y + REAL(xp) :: pade_numerator_coeffs(N_NUM + 1), pade_denominator_coeffs(N_NUM + 1) + REAL(xp) :: numerator_sum + REAL(xp) :: denominator_sum + INTEGER :: m + ! If N_NUM == 0, then the approximation should be the same as the Taylor approx. of N_NUM: + IF (N_NUM == 0) THEN + pade_kernel_n = taylor_kernel_n(N_NUM, n, y) + ELSE + CALL find_pade_coefficients(pade_numerator_coeffs, pade_denominator_coeffs, n, N_NUM, N_DEN) + numerator_sum = 0 + denominator_sum = 0 + DO m = 0, N_NUM + numerator_sum = numerator_sum + pade_numerator_coeffs(m + 1) * y ** m + END DO + DO m = 0, N_NUM + denominator_sum = denominator_sum + pade_denominator_coeffs(m + 1) * y ** m + END DO + pade_kernel_n = numerator_sum / denominator_sum + END IF +END FUNCTION pade_kernel_n + + +SUBROUTINE find_pade_coefficients(pade_num_coeffs, pade_denom_coeffs, n, N_NUM, N_DEN) + IMPLICIT NONE +#ifdef SINGLE_PRECISION + EXTERNAL :: SGESV ! Use DGESV rather than SGESV for double precision +#else + EXTERNAL :: DGESV ! Use DGESV rather than SGESV for double precision +#endif + INTEGER, INTENT (IN) :: n, N_NUM, N_DEN ! index of the considered Kernel + REAL(xp), INTENT(OUT) :: pade_num_coeffs(N_NUM + 1), pade_denom_coeffs(N_NUM + 1) ! OUT rather than INOUT to make sure no information is retained from previous Kernel computations + + REAL(xp) :: taylor_kernel_coeffs(N_NUM + N_NUM + 1), denom_matrix(N_NUM, N_NUM) + INTEGER :: m, j + REAL(xp) :: m_dp, n_dp + INTEGER :: return_code ! for DGESV + REAL(xp) :: pivot(N_NUM) ! for DGESV + + n_dp = REAL(n, xp) + + ! First find the kernel Taylor expansion coefficients + DO m = 0, N_NUM + N_NUM ! m here counts the order of the derivatives + m_dp = REAL(m, xp) + + IF (m < n) THEN + taylor_kernel_coeffs(m + 1) = 0 + ELSE + taylor_kernel_coeffs(m + 1) = (-1)**(n + m) / (GAMMA(n_dp + 1._xp) * GAMMA(m_dp - n_dp + 1._xp)) + END IF + END DO + + ! Next construct the denominator solving matrix + DO m = 1, N_NUM + DO j = 1, N_NUM + IF (N_NUM + m - j < 0) THEN + denom_matrix(m, j) = 0 + ELSE + denom_matrix(m, j) = taylor_kernel_coeffs(N_NUM + m - j + 1) + END IF + END DO + END DO + + ! Then solve for the denominator coefficients, setting the first one to 1 + !!!! SOLVER NOT YET IMPLEMENTED!!!! + pade_denom_coeffs(1) = 1 + pade_denom_coeffs(2:) = - taylor_kernel_coeffs(N_NUM + 2 : N_NUM + N_NUM + 1) ! First acts as RHS vector for equation, is then transformed to solution by DGESV +#ifdef SINGLE_PRECISION + CALL SGESV(N_NUM, 1, denom_matrix, N_NUM, pivot, pade_denom_coeffs(2:), N_NUM, return_code) ! LAPACK solver for matrix equation. Note that denom_matrix is now no longer as expected +#else + CALL DGESV(N_NUM, 1, denom_matrix, N_NUM, pivot, pade_denom_coeffs(2:), N_NUM, return_code) ! LAPACK solver for matrix equation. Note that denom_matrix is now no longer as expected +#endif + + + ! Print an error message in case there was a problem with the solver + IF (return_code /= 0) THEN + WRITE (*,*) 'An error occurred in the solving for the Pade denominator coefficients. The error code is: ', return_code + END IF + + ! Finally compute the numerator coefficients + DO m = 1, N_NUM + 1 + pade_num_coeffs(m) = 0 ! As the array is not automatically filled with zeros + DO j = 1, m + !num_matrix(m, j) = taylor_kernel_coeffs(m - j + 1) + pade_num_coeffs(m) = pade_num_coeffs(m) + pade_denom_coeffs(j) * taylor_kernel_coeffs(m - j + 1) + END DO + END DO + +END SUBROUTINE find_pade_coefficients +!******************************************************************************! + + END MODULE numerics diff --git a/src/solve_EM_fields.F90 b/src/solve_EM_fields.F90 index d28a4aca89bd83c3e024baa31cc9e580e563d9a6..dfbc9d2dedfa97ca32b572316832535855cf0fec 100644 --- a/src/solve_EM_fields.F90 +++ b/src/solve_EM_fields.F90 @@ -69,10 +69,7 @@ CONTAINS !!!!!!!!!!!!!!! adiabatic ion model ! Candy et al. 2007, rho_i = -q_i/tau_i phi IF (ADIAB_I) THEN - DO iz = 1,local_nz - izi = iz+ngz/2 - rho(iz) = rho(iz) - 0*q_i/tau_i * phi(iky,ikx,izi) - ENDDO + ! No ion contribution ENDIF !!!!!!!!!!!!!!! Inverting the poisson equation diff --git a/src/species_mod.F90 b/src/species_mod.F90 index 2d86e0601d40fb81ea8843927a1757a225afdb0c..43c8557bbb2176aaccd4874f8a0de76b11e6834b 100644 --- a/src/species_mod.F90 +++ b/src/species_mod.F90 @@ -10,6 +10,7 @@ MODULE species REAL(xp) :: q_ ! Charge REAL(xp) :: k_N_ ! density drive (L_ref/L_Ni) REAL(xp) :: k_T_ ! temperature drive (L_ref/L_Ti) + REAL(xp) :: mu_ ! species-wise hyperdiffusion (not tested) !! Arrays to store all species features CHARACTER(len=32),& ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: name ! name of the species @@ -19,6 +20,7 @@ MODULE species REAL(xp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: k_N ! density drive (L_ref/L_Ni) REAL(xp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: k_T ! temperature drive (L_ref/L_Ti) REAL(xp), ALLOCATABLE, DIMENSION(:,:),PUBLIC, PROTECTED :: nu_ab ! Collision frequency tensor + REAL(xp), ALLOCATABLE, DIMENSION(:), PUBLIC, PROTECTED :: mu ! Hyperdiffusion !! Auxiliary variables to store precomputation REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: tau_q ! factor of the magnetic moment coupling REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: q_tau ! charge/temp ratio @@ -43,7 +45,7 @@ CONTAINS INTEGER :: ia,ib ! expected namelist in the input file NAMELIST /SPECIES/ & - name_, tau_, sigma_, q_, k_N_, k_T_ + name_, tau_, sigma_, q_, k_N_, k_T_, mu_ ! allocate the arrays of species parameters CALL species_allocate ! loop over the species namelists in the input file @@ -56,6 +58,7 @@ CONTAINS q_ = 1._xp k_N_ = 2.22_xp k_T_ = 6.96_xp + mu_ = 0._xp ! read input READ(lu_in,species) ! place values found in the arrays @@ -66,6 +69,7 @@ CONTAINS k_N(ia) = k_N_ k_T(ia) = k_T_ tau_q(ia) = tau_/q_ + mu(ia) = mu_ ! precompute factors q_tau(ia) = q_/tau_ sqrtTau_q(ia) = sqrt(tau_)/q_ @@ -130,6 +134,7 @@ CONTAINS CALL attach(fid, TRIM(str), "q", q(ia)) CALL attach(fid, TRIM(str), "k_N", k_N(ia)) CALL attach(fid, TRIM(str), "k_T", k_T(ia)) + CALL attach(fid, TRIM(str), "mu", mu(ia)) ENDDO END SUBROUTINE species_outputinputs @@ -144,6 +149,7 @@ CONTAINS ALLOCATE( q(Na)) ALLOCATE( k_N(Na)) ALLOCATE( k_T(Na)) + ALLOCATE( mu(Na)) ALLOCATE( tau_q(Na)) ALLOCATE( q_tau(Na)) ALLOCATE( sqrtTau_q(Na)) diff --git a/testcases/ITG_zpinch/3D/fort_00.90 b/testcases/ITG_zpinch/3D/fort_00.90 new file mode 100644 index 0000000000000000000000000000000000000000..d6255d0b0adb301845994752276cbcb0c4646089 --- /dev/null +++ b/testcases/ITG_zpinch/3D/fort_00.90 @@ -0,0 +1,104 @@ +&BASIC + nrun = 99999999 + dt = 0.01 + tmax = 500 + maxruntime = 72000 + job2load = -1 +/ +&GRID + pmax = 2 + jmax = 1 + Nx = 128 + Lx = 80 + Ny = 64 + Ly = 80 + Nz = 24 + SG = .t. + Nexc = 1 +/ +&GEOMETRY + geom = 'Z-pinch' + q0 = 0.0 + shear = 0.0 + eps = 0.0 + kappa = 1.0 + s_kappa= 0.0 + delta = 0.0 + s_delta= 0.0 + zeta = 0.0 + s_zeta = 0.0 + parallel_bc = 'shearless' + shift_y= 0.0 + Npol = 1 +/ +&DIAGNOSTICS + dtsave_0d = 1 + dtsave_1d = -1 + dtsave_2d = -1 + dtsave_3d = 0.5 + dtsave_5d = 50 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'nonlinear' + Na = 1 ! number of species + mu_x = 1.0 + mu_y = 1.0 + N_HD = 4 + mu_z = 1.0 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.05 + beta = 0.0 + ADIAB_E = .t. +/ +&CLOSURE + !hierarchy_closure='truncation' + hierarchy_closure='max_degree' + dmax = 2 + nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing) + nmax = 1 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 0.0 + k_T_ = 6.96 +/ +&SPECIES + ! electrons + name_ = 'electrons' + tau_ = 1.0 + sigma_= 0.023338 + q_ =-1.0 + k_N_ = 1.6 + k_T_ = 0.4 +/ +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .f. + INTERSPECIES = .true. + mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5' +/ +&INITIAL + INIT_OPT = 'blob' !(phi,blob) + ACT_ON_MODES = 'donothing' + init_background = 0.0 + init_noiselvl = 0.005 + iseed = 42 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' +/ diff --git a/testcases/ITG_zpinch/NO_SG/fort_00.90 b/testcases/ITG_zpinch/NO_SG/fort_00.90 new file mode 100644 index 0000000000000000000000000000000000000000..381ace207c4e3c115a4066627e74a2ad529a10fb --- /dev/null +++ b/testcases/ITG_zpinch/NO_SG/fort_00.90 @@ -0,0 +1,103 @@ +&BASIC + nrun = 99999999 + dt = 0.02 + tmax = 100 + maxruntime = 72000 + job2load = -1 +/ +&GRID + pmax = 2 + jmax = 1 + Nx = 64 + Lx = 80 + Ny = 48 + Ly = 80 + Nz = 16 + SG = .f. + Nexc = 1 +/ +&GEOMETRY + geom = 'Z-pinch' + q0 = 0.0 + shear = 0.0 + eps = 0.0 + kappa = 1.0 + s_kappa= 0.0 + delta = 0.0 + s_delta= 0.0 + zeta = 0.0 + s_zeta = 0.0 + parallel_bc = 'shearless' + shift_y= 0.0 +/ +&DIAGNOSTICS + dtsave_0d = 1 + dtsave_1d = -1 + dtsave_2d = -1 + dtsave_3d = 0.5 + dtsave_5d = 50 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'nonlinear' + Na = 1 ! number of species + mu_x = 0.0 + mu_y = 0.0 + N_HD = 4 + mu_z = 0.0 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.05 + beta = 0.0 + ADIAB_E = .t. +/ +&CLOSURE + !hierarchy_closure='truncation' + hierarchy_closure='max_degree' + dmax = 2 + nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing) + nmax = -1 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 0.0 + k_T_ = 4.0 +/ +&SPECIES + ! electrons + name_ = 'electrons' + tau_ = 1.0 + sigma_= 0.023338 + q_ =-1.0 + k_N_ = 1.6 + k_T_ = 0.4 +/ +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .t. + INTERSPECIES = .true. + mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5' +/ +&INITIAL + INIT_OPT = 'phi' !(phi,blob) + ACT_ON_MODES = 'donothing' + init_background = 0.0 + init_noiselvl = 0.005 + iseed = 42 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' +/ diff --git a/testcases/ITG_zpinch/PAR_DIFF/fort_00.90 b/testcases/ITG_zpinch/PAR_DIFF/fort_00.90 new file mode 100644 index 0000000000000000000000000000000000000000..30fdace3f704b6418134cc0d8df9c2992bef1716 --- /dev/null +++ b/testcases/ITG_zpinch/PAR_DIFF/fort_00.90 @@ -0,0 +1,103 @@ +&BASIC + nrun = 99999999 + dt = 0.02 + tmax = 100 + maxruntime = 72000 + job2load = -1 +/ +&GRID + pmax = 2 + jmax = 1 + Nx = 64 + Lx = 80 + Ny = 48 + Ly = 80 + Nz = 16 + SG = .f. + Nexc = 1 +/ +&GEOMETRY + geom = 'Z-pinch' + q0 = 0.0 + shear = 0.0 + eps = 0.0 + kappa = 1.0 + s_kappa= 0.0 + delta = 0.0 + s_delta= 0.0 + zeta = 0.0 + s_zeta = 0.0 + parallel_bc = 'shearless' + shift_y= 0.0 +/ +&DIAGNOSTICS + dtsave_0d = 1 + dtsave_1d = -1 + dtsave_2d = -1 + dtsave_3d = 0.5 + dtsave_5d = 50 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'nonlinear' + Na = 1 ! number of species + mu_x = 0.0 + mu_y = 0.0 + N_HD = 4 + mu_z = 0.2 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.05 + beta = 0.0 + ADIAB_E = .t. +/ +&CLOSURE + !hierarchy_closure='truncation' + hierarchy_closure='max_degree' + dmax = 2 + nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing) + nmax = -1 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 0.0 + k_T_ = 4.0 +/ +&SPECIES + ! electrons + name_ = 'electrons' + tau_ = 1.0 + sigma_= 0.023338 + q_ =-1.0 + k_N_ = 1.6 + k_T_ = 0.4 +/ +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .t. + INTERSPECIES = .true. + mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5' +/ +&INITIAL + INIT_OPT = 'phi' !(phi,blob) + ACT_ON_MODES = 'donothing' + init_background = 0.0 + init_noiselvl = 0.005 + iseed = 42 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' +/ diff --git a/testcases/ITG_zpinch/fort.90 b/testcases/ITG_zpinch/fort.90 deleted file mode 100644 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000 diff --git a/testcases/ITG_zpinch/fort_00.90 b/testcases/ITG_zpinch/fort_00.90 index 87253b73dd8669e3481ffc76c1e920fd59f5d150..522164900fc4152cfd9b103684cec6ee97433096 100644 --- a/testcases/ITG_zpinch/fort_00.90 +++ b/testcases/ITG_zpinch/fort_00.90 @@ -1,19 +1,19 @@ &BASIC nrun = 99999999 dt = 0.01 - tmax = 250 + tmax = 500 maxruntime = 72000 job2load = -1 / &GRID pmax = 2 jmax = 1 - Nx = 64 - Lx = 120 - Ny = 48 - Ly = 120 - Nz = 16 - SG = .f. + Nx = 128 + Lx = 80 + Ny = 64 + Ly = 80 + Nz = 1 + SG = .t. Nexc = 1 / &GEOMETRY @@ -34,8 +34,8 @@ dtsave_0d = 1 dtsave_1d = -1 dtsave_2d = -1 - dtsave_3d = 1 - dtsave_5d = 10 + dtsave_3d = 0.5 + dtsave_5d = 50 write_doubleprecision = .f. write_gamma = .t. write_hf = .t. @@ -49,8 +49,8 @@ &MODEL LINEARITY = 'nonlinear' Na = 1 ! number of species - mu_x = 0.0 - mu_y = 0.0 + mu_x = 1.0 + mu_y = 1.0 N_HD = 4 mu_z = 0.0 HYP_V = 'hypcoll' @@ -59,14 +59,13 @@ nu = 0.05 beta = 0.0 ADIAB_E = .t. - tau_e = 1.0 / &CLOSURE !hierarchy_closure='truncation' hierarchy_closure='max_degree' dmax = 2 nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing) - nmax = -1 + nmax = 1 / &SPECIES ! ions @@ -75,7 +74,7 @@ sigma_= 1.0 q_ = 1.0 k_N_ = 0.0 - k_T_ = 4.0 + k_T_ = 6.96 / &SPECIES ! electrons @@ -88,7 +87,7 @@ / &COLLISION collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) - GK_CO = .t. + GK_CO = .f. INTERSPECIES = .true. mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5' / diff --git a/testcases/ITG_zpinch/fort_01.90 b/testcases/ITG_zpinch/fort_01.90 new file mode 100644 index 0000000000000000000000000000000000000000..3f22e004be173a628ab9e2e72f95ef8eda9446d1 --- /dev/null +++ b/testcases/ITG_zpinch/fort_01.90 @@ -0,0 +1,103 @@ +&BASIC + nrun = 99999999 + dt = 0.01 + tmax = 1000 + maxruntime = 72000 + job2load = 0 +/ +&GRID + pmax = 2 + jmax = 1 + Nx = 128 + Lx = 80 + Ny = 64 + Ly = 80 + Nz = 16 + SG = .t. + Nexc = 1 +/ +&GEOMETRY + geom = 'Z-pinch' + q0 = 0.0 + shear = 0.0 + eps = 0.0 + kappa = 1.0 + s_kappa= 0.0 + delta = 0.0 + s_delta= 0.0 + zeta = 0.0 + s_zeta = 0.0 + parallel_bc = 'shearless' + shift_y= 0.0 +/ +&DIAGNOSTICS + dtsave_0d = 1 + dtsave_1d = -1 + dtsave_2d = -1 + dtsave_3d = 0.5 + dtsave_5d = 50 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'nonlinear' + Na = 1 ! number of species + mu_x = 1.0 + mu_y = 1.0 + N_HD = 4 + mu_z = 0.0 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.05 + beta = 0.0 + ADIAB_E = .t. +/ +&CLOSURE + !hierarchy_closure='truncation' + hierarchy_closure='max_degree' + dmax = 2 + nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing) + nmax = 1 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 0.0 + k_T_ = 6.96 +/ +&SPECIES + ! electrons + name_ = 'electrons' + tau_ = 1.0 + sigma_= 0.023338 + q_ =-1.0 + k_N_ = 1.6 + k_T_ = 0.4 +/ +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .f. + INTERSPECIES = .true. + mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5' +/ +&INITIAL + INIT_OPT = 'blob' !(phi,blob) + ACT_ON_MODES = 'donothing' + init_background = 0.0 + init_noiselvl = 0.005 + iseed = 42 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' +/ diff --git a/testcases/Rosenbluth_Hinton_test/fort_00.90 b/testcases/Rosenbluth_Hinton_test/fort_00.90 new file mode 100644 index 0000000000000000000000000000000000000000..bc86583eef35435c1298ac6ad684f9aaa7c8bfed --- /dev/null +++ b/testcases/Rosenbluth_Hinton_test/fort_00.90 @@ -0,0 +1,101 @@ +&BASIC + nrun = 99999999 + dt = 0.025 + tmax = 25 + maxruntime = 72000 + job2load = -1 +/ +&GRID + pmax = 32 + jmax = 0 + Nx = 2 + Lx = 120 + Ny = 2 + Ly = 120 + Nz = 24 + SG = .f. + Nexc = 0 +/ +&GEOMETRY + geom = 's-alpha' + !geom = 'miller' + q0 = 1.4 + shear = 0.0 + eps = 0.1 + kappa = 1.0 + s_kappa= 0.0 + delta = 0.0 + s_delta= 0.0 + zeta = 0.0 + s_zeta = 0.0 + parallel_bc = 'dirichlet' + shift_y= 0.0 +/ +&DIAGNOSTICS + dtsave_0d = 1 + dtsave_1d = -1 + dtsave_2d = -1 + dtsave_3d = 0.1 + dtsave_5d = 20 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'linear' + Na = 1 ! number of species + mu_x = 0.0 + mu_y = 0.0 + N_HD = 4 + mu_z = 0.0 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.01 + beta = 0.0 + ADIAB_E = .t. +/ +&CLOSURE + hierarchy_closure='truncation' + dmax = -1 + nonlinear_closure='truncation' + nmax = 0 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 0 + k_T_ = 0 +/ + +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .f. + INTERSPECIES = .true. + mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5' +/ +&INITIAL + INIT_OPT = 'mom00' + !INIT_OPT = 'mom00_mode' + init_background = 1.0 + init_noiselvl = 0 + iseed = 42 + Nmodes = 1 +/ +&MODE + I_ = 1 + J_ = 0 + amp_ = 1.0 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' +/ diff --git a/testcases/cyclone_example/Z_pinch/fort_00.90 b/testcases/cyclone_example/Z_pinch/fort_00.90 new file mode 100644 index 0000000000000000000000000000000000000000..96c03136dd6b2778c54337686dfa9457e0050baf --- /dev/null +++ b/testcases/cyclone_example/Z_pinch/fort_00.90 @@ -0,0 +1,96 @@ +&BASIC + nrun = 99999999 + dt = 0.01 + tmax = 50 + maxruntime = 72000 + job2load = -1 +/ +&GRID + pmax = 4 + jmax = 1 + Nx = 64 + Lx = 120 + Ny = 48 + Ly = 120 + Nz = 16 + SG = .f. + Nexc = 0 +/ +&GEOMETRY + geom = 'Z-pinch' + !geom = 'miller' + q0 = 100!1.4 + shear = 0.0!0.8 + eps = 0.001!0.18 + kappa = 1.0 + s_kappa= 0.0 + delta = 0.0 + s_delta= 0.0 + zeta = 0.0 + s_zeta = 0.0 + parallel_bc = 'dirichlet' + shift_y= 0.0 +/ +&DIAGNOSTICS + dtsave_0d = 0.1 + dtsave_1d = -1 + dtsave_2d = -1 + dtsave_3d = 1 + dtsave_5d = 20 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'nonlinear' + Na = 1 ! number of species + mu_x = 0.0 + mu_y = 0.0 + N_HD = 4 + mu_z = 1.0 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.05 + beta = 0.0 + ADIAB_E = .t. +/ +&CLOSURE + hierarchy_closure='truncation' + dmax = -1 + nonlinear_closure='truncation' + nmax = 0 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 2.22 + k_T_ = 6.96 +/ + +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .f. + INTERSPECIES = .true. + mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5' +/ +&INITIAL + INIT_OPT = 'blob' + ACT_ON_MODES = 'donothing' + init_background = 0.0 + init_noiselvl = 0.005 + iseed = 42 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' + !numerical_scheme = 'SSP_RK3' +/ diff --git a/testcases/cyclone_example/Z_pinch_limit/fort_00.90 b/testcases/cyclone_example/Z_pinch_limit/fort_00.90 new file mode 100644 index 0000000000000000000000000000000000000000..a555316f0a3503c4e6ec430e7cb5ca90c4fad25e --- /dev/null +++ b/testcases/cyclone_example/Z_pinch_limit/fort_00.90 @@ -0,0 +1,96 @@ +&BASIC + nrun = 99999999 + dt = 0.01 + tmax = 0.1 + maxruntime = 72000 + job2load = -1 +/ +&GRID + pmax = 4 + jmax = 1 + Nx = 64 + Lx = 120 + Ny = 48 + Ly = 120 + Nz = 16 + SG = .f. + Nexc = 0 +/ +&GEOMETRY + !geom = 's-alpha' + geom = 'miller' + q0 = 1000!1.4 + shear = 0.0!0.8 + eps = 0.001!0.18 + kappa = 1.0 + s_kappa= 0.0 + delta = 0.0 + s_delta= 0.0 + zeta = 0.0 + s_zeta = 0.0 + parallel_bc = 'dirichlet' + shift_y= 0.0 +/ +&DIAGNOSTICS + dtsave_0d = 0.1 + dtsave_1d = -1 + dtsave_2d = -1 + dtsave_3d = 1 + dtsave_5d = 20 + write_doubleprecision = .f. + write_gamma = .t. + write_hf = .t. + write_phi = .t. + write_Na00 = .t. + write_Napj = .t. + write_dens = .t. + write_fvel = .t. + write_temp = .t. +/ +&MODEL + LINEARITY = 'nonlinear' + Na = 1 ! number of species + mu_x = 0.0 + mu_y = 0.0 + N_HD = 4 + mu_z = 1.0 + HYP_V = 'hypcoll' + mu_p = 0.0 + mu_j = 0.0 + nu = 0.05 + beta = 0.0 + ADIAB_E = .t. +/ +&CLOSURE + hierarchy_closure='truncation' + dmax = -1 + nonlinear_closure='truncation' + nmax = 0 +/ +&SPECIES + ! ions + name_ = 'ions' + tau_ = 1.0 + sigma_= 1.0 + q_ = 1.0 + k_N_ = 2.22 + k_T_ = 6.96 +/ + +&COLLISION + collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) + GK_CO = .f. + INTERSPECIES = .true. + mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5' +/ +&INITIAL + INIT_OPT = 'blob' + ACT_ON_MODES = 'donothing' + init_background = 0.0 + init_noiselvl = 0.005 + iseed = 42 +/ +&TIME_INTEGRATION + numerical_scheme = 'RK4' + !numerical_scheme = 'SSP_RK3' +/ diff --git a/testcases/cyclone_example/fort_00.90 b/testcases/cyclone_example/fort_00.90 index ce02b619f0d1394c95dc6ecde45bdf99934531b2..abeb5052a73822c2351f43a623f63c2bd25f27b6 100644 --- a/testcases/cyclone_example/fort_00.90 +++ b/testcases/cyclone_example/fort_00.90 @@ -1,7 +1,7 @@ &BASIC nrun = 99999999 dt = 0.01 - tmax = 5 + tmax = 50 maxruntime = 72000 job2load = -1 / @@ -53,7 +53,7 @@ mu_x = 0.0 mu_y = 0.0 N_HD = 4 - mu_z = 1.0 + mu_z = 2.0 HYP_V = 'hypcoll' mu_p = 0.0 mu_j = 0.0 @@ -79,7 +79,7 @@ &COLLISION collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau) - GK_CO = .f. + GK_CO = .t. INTERSPECIES = .true. mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5' / diff --git a/wk/Ch7_HF_analysis.m b/wk/Ch7_HF_analysis.m new file mode 100644 index 0000000000000000000000000000000000000000..750ec4b91f6f3a950b95b72ff7f7e2d61ba75ca1 --- /dev/null +++ b/wk/Ch7_HF_analysis.m @@ -0,0 +1,45 @@ +PARTITION = '/misc/gyacomo23_outputs/paper_3/'; +switch 3 +case 1 + SIM_SET_NAME = 'Multi-scale'; + E_FLUX = 1; + FOLDER = 'DIIID_fullphys_rho95_geom_scan/multi_scale/multi_scale_3x2x768x192x24'; + % FOLDER = 'DIIID_fullphys_rho95_geom_scan/multi_scale/multi_scale_5x2x768x192x24'; +case 2 + SIM_SET_NAME = 'Ion-scale'; + E_FLUX = 1; + FOLDER = 'DIIID_fullphys_rho95_geom_scan/ion_scale/ion_scale_5x2x256x64x32_tau_1_RN_0'; +case 3 + SIM_SET_NAME = 'Adiab. e.'; + E_FLUX = 0; + FOLDER = 'DIIID_adiab_e_rho95_geom_scan/5x2x256x64x32_tau_1_RN_0/'; +case 4 + SIM_SET_NAME = 'Adiab. e.'; + E_FLUX = 0; + FOLDER = 'DIIID_adiab_e_rho95_geom_scan/5x2x256x64x32_tau_1_RN_0/'; +end + +GEOM = {'NT','0T','PT'}; +COLO = {'b','k','r'}; +figure +for i = 1:3 + OPTION = GEOM{i}; + datadir = [PARTITION,FOLDER,'/',OPTION]; + out = read_flux_out_XX(datadir,0); + + [ts, is] = sort(out.t); + Pxis = out.Pxi(is); + Qxis = out.Qxi(is); + + if E_FLUX + Pxes = out.Pxe(is); + Qxes = out.Qxe(is); + subplot(211) + plot(ts,Qxes,COLO{i}); hold on; + subplot(212) + title(SIM_SET_NAME) + end + plot(ts,Qxis,COLO{i}); hold on; +end +title(SIM_SET_NAME) + diff --git a/wk/HEL_lin_solver.m b/wk/HEL_GM_eig_solver.m similarity index 82% rename from wk/HEL_lin_solver.m rename to wk/HEL_GM_eig_solver.m index 752fb87f5e8c45cd9aff50f6f19cc2df5f8a8f92..457f4417847a408aead1211f6273028cd603cc04 100644 --- a/wk/HEL_lin_solver.m +++ b/wk/HEL_GM_eig_solver.m @@ -1,24 +1,33 @@ % We formulate the linear system as dN/dt + A*N = 0 +% parameters +q = 1; % Safety factor +Jxyz = 1; % Jacobian +hatB = 1; % normalized B field +dzlnB = 0; % B parallel derivative +ddz = 0; % parallel derivative operator +tau = 0.001; % ion-electron temperature ratio +IVAN = 0; % flag for Ivanov scaling +kT = 7; % scaled temperature gradient +chi = 0.16; % scaled collisionality +kx_a = 0;linspace(-2.0,2.0,256); % radial fourier grid +ky_a = linspace(0.01,5,256); % binormal fourier grid + +% translate into parameters in GM formulation +RN = 0; +RT = kT*2*q/tau; +NU = chi*3/2/tau/(4-IVAN*2.25); +MU = 0.0; -q = 1; -Jxyz = 1; hatB = 1; dzlnB = 0; -ddz = 0; -tau = 0.001; -IVAN = 1; %ordering O1_n = 1; O1_u = 1-IVAN; O1_Tpar = 1-IVAN; O1_Tper = 1-IVAN; -RN = 0; -RT = 1*2*q/tau; -NU = 0.1*3/2/tau/(4-IVAN*2.25); -MU = 0.0; - -kx_a = 0;linspace(-2.0,2.0,256); -ky_a = linspace(0.01,3.5,256); +% array of most unstable growth rates and corresponding frequencies g_a = zeros(numel(kx_a),numel(ky_a)); w_a = g_a; + +% Evaluation of the matrices and solver for i = 1:numel(kx_a) for j = 1:numel(ky_a) kx = kx_a(i); @@ -116,26 +125,8 @@ for j = 1:numel(ky_a) CDG(4,3) =-2/3*sqrt(2)*K0*(K0-2*K1*O1_Tper); CDG(4,4) = 4/3*(K0-2*K1*O1_Tper)^2 + 2*tau*(K0-K1)^2*lperp*O1_Tper; CDG(4,5) =-2/3*q/tau*(K0-K1)*(3*tau*K0^2*lperp-K0*K1*(4+3*tau*lperp)+K1^2*(8+3*tau*lperp)); - if 1 - C = NU*(CLB + CDG); - else %to check - C = zeros(4,5); - C(1,1) = -2/3*tau*lperp*(4*tau*lperp); - C(1,3) = -2/3*tau*lperp*(sqrt(2)-2*sqrt(2)*tau*lperp); - C(1,4) = -2/3*tau*lperp*(1-tau*lperp); - C(1,5) = -2/3*tau*lperp*(5*q*lperp - 11*q*tau*lperp^2); - C(2,2) = -(1+2*tau*lperp); - C(3,1) = -2/3*(sqrt(2)*tau*lperp); - C(3,3) = -2/3*(2+5*tau*lperp); - C(3,4) = -2/3*(sqrt(2)-4*sqrt(2)*tau*lperp); - C(3,5) = -2/3*(2*sqrt(2)*q*lperp+8*sqrt(2)*q*tau*lperp^2); - C(4,1) = -2/3*(tau*lperp - tau^2*lperp^2); - C(4,3) = -2/3*(sqrt(2)*4*sqrt(2)*tau*lperp); - C(4,4) = -2/3*(1+12*tau*lperp); - C(4,5) = -2/3*(2*q*lperp+9*q*tau*lperp^2); - C = NU*C; + C = NU*(CLB + CDG); - end % Numerical diffusion Diff = MU*lperp^2*eye(4); @@ -166,7 +157,11 @@ for j = 1:numel(ky_a) gamma =-real(lambda); omega = imag(lambda); [g_a(i,j),idx] = max(gamma); - % w_a(i,j) = omega(idx); + % gsorted = sort(gamma); + % g_a(i,j) = gsorted(1); + w_a(i,j) = omega(idx); + % wsorted = sort(omega); + % g_a(i,j) = wsorted(1); w_a(i,j) = max(omega); end end diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m index 215a920d5753933fa524a6cdb26377480ee81713..4e0e5814654e347f9d832e3af220307ef7567f15 100644 --- a/wk/fast_analysis.m +++ b/wk/fast_analysis.m @@ -7,7 +7,7 @@ default_plots_options % Partition of the computer where the data have to be searched % PARTITION='/Users/ahoffmann/gyacomo/results/paper_3/'; PARTITION='/misc/gyacomo23_outputs/paper_3/'; -% PARTITION = ''; +% PARTITION = '../results/paper_3/'; %% Paper 3 % resdir = 'DTT_rho85/3x2x192x48x32'; % resdir = 'DTT_rho85/3x2x192x48x32_NT'; @@ -17,88 +17,168 @@ PARTITION='/misc/gyacomo23_outputs/paper_3/'; % resdir = 'LM_DIIID_rho95/3x2x512x92x32'; % resdir = 'DIIID_LM_rho90/3x2x256x128x32'; % resdir = 'DTT_rho85_geom_scan/P8_J4_delta_nuDGGK_conv_test/delta_-0.3_nu_0.9'; -resdir = 'NT_DIIID_Austin2019_rho95/3x2x256x64x32'; -% resdir = '../testcases/cyclone_example'; +% resdir = 'NT_DIIID_Austin2019_rho95/3x2x256x64x32'; + +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/largerbox_moremodes'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/PT/3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/NT/3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/VNT/3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/VPT/3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/CIRC/3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/ELONG/3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/PT/lin_3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/NT/lin_3x2x128x32x32'; + +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/huge'; + +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_1000/PT/3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_1000/NT/3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_1000/PT/lin_3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_1000/NT/lin_3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_750/PT/lin_3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_750/NT/lin_3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_250/PT/lin_3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_250/NT/lin_3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-3/huge'; + +% resdir = 'DIIID_rho95_cold_ions_tau1e-2/PT/3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e-2/NT/3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e0/PT/3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e0/PT/5x3x192x48x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e0/NT/3x2x128x32x32'; +% resdir = 'DIIID_rho95_cold_ions_tau1e0/NT/5x3x192x48x32'; +% resdir = 'DTT_rho85_cold_ions_tau1e-3/NT/3x2x128x32x32'; +% % resdir = '../testcases/cyclone_example'; % resdir = '../testcases/CBC_ExBshear/std'; % resdir = '../results/paper_3/HM_DTT_rho98/3x2x128x64x64'; - %% -J0 = 00; J1 = 10; -% Load basic info (grids and time traces) +% resdir = 'DIIID_cold_ions_rho95_geom_scan/3x2x192x48x32_RT_1000_eps_q0_scan/NT/eps_0.35_q0_4.0'; +% resdir = 'DTT_rho85_geom_scan/P2_J1_PT_sfact_shear_scan/shear_2.7_q0_-2.9'; +% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_kT_3500/128x32x24'; +% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_kT_3500/256x64x24'; +% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_colless'; +% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_kT_2000/lin_128x32x24'; +% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_kT_2000/128x32x24'; +% PARTITION = ''; resdir ='../results/HEL_CBC/CBC_21/128x32x24'; +% PARTITION = ''; resdir ='../results/HEL_CBC/192x48x24'; +% PARTITION = ''; resdir ='../testcases/Ivanov_2020'; +% PARTITION = ''; resdir ='../testcases/Hasegawa_Wakatani'; +% resdir = 'HEL_CBC/256x92x24_max_trunc'; + +% resdir ='HEL_CBC/192x48x24'; +% resdir ='HEL_CBC/256x92x24'3; +% resdir ='HEL_CBC/256x256x32/k_N__0.0_k_T__1750'; + +PARTITION = '/misc/gyacomo23_outputs/paper_3/DIIID_rho_95_Tstudy/'; +% resdir = 'multi_scale_3x2x512x128x24'; +% resdir = 'multi_scale_3x2x512x128x24_larger_box'; +% resdir = 'multi_scale_3x2x768x192x24/continue_with_gradN_and_tau'; +% resdir = 'multi_scale/multi_scale_5x2x768x192x24/PT'; +% resdir = 'NT'; +% resdir = 'electron_scale_3x2x256x64x24'; +% resdir = 'electron_scale_3x2x128x64x24'; +% resdir = 'ion_scale_3x2x192x48x32_larger_box'; +% resdir = 'ion_scale_3x2x256x64x32'; % PT and NT also +% resdir = 'ion_scale/ion_scale_5x2x256x64x32_tau_1_RN_0/NT'; +% resdir = 'adiab_e/5x2x256x64x32_tau_1_RN_0/NT'; + +resdir = 'multi_scale/multi_scale_3x2x768x192x24/NT'; +% resdir = 'multi_scale/multi_scale_3x2x512x128x24_larger_box/NT'; +% resdir = 'multi_scale/multi_scale_3x2x768x192x24/continue_with_gradN_and_tau'; +% resdir = 'ion_scale/ion_scale_5x2x256x64x32_tau_1_RN_0/NT'; +% resdir = 'adiab_e/5x2x256x64x32_tau_1_RN_0/0T'; +% resdir = 'hot_electrons/hot_electrons_256x64x24/0T'; + + + +% PARTITION = '/misc/gyacomo23_outputs/paper_3/'; +% resdir = 'DIIID_HEL_rho95/PT'; + DATADIR = [PARTITION,resdir,'/']; +%% +J0 = 00; J1 = 20; + +% Load basic info (grids and time traces) data = {}; data = compile_results_low_mem(data,DATADIR,J0,J1); - +[data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00'); +data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +try +data.Ne00 = reshape(data.Na00(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +catch +end +[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi'); +if 1 + %% +[data.TEMP, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'temp'); +% [data.UPAR, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'upar'); +[data.DENS, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'dens'); +data.TEMP_I = reshape(data.TEMP(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +% data.UPAR_I = reshape(data.UPAR(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +data.DENS_I = reshape(data.DENS(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +if data.inputs.Na > 1 + data.TEMP_E = reshape(data.TEMP(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); + data.DENS_E = reshape(data.DENS(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); + data.Ne00 = reshape(data.Na00(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +end +end if 1 %% Plot transport and phi radial profile -[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi'); +% [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi'); % [data.PSI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'psi'); -options.TAVG_0 = 100; -options.TAVG_1 = 1000; +options.TAVG_0 = data.Ts3D(end)/2; +options.TAVG_1 = data.Ts3D(end); options.NCUT = 5; % Number of cuts for averaging and error estimation 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 = 'n_e'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) +% options.ST_FIELD = 'u_i'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) +% options.ST_FIELD = 'T_i'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) +% options.ST_FIELD = 'n_i T_i'; % 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.ST_FIELD = 'Q_{xi}'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) +% options.ST_FIELD = 'G_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) +% options.ST_FIELD = 'w_{Ez}'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) +% options.ST_FIELD = 'v_{Ey}'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) +% options.ST_FIELD = 'N_i^{00}'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) options.INTERP = 0; options.RESOLUTION = 256; plot_radial_transport_and_spacetime(data,options); end -if 0 -%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Options -[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi'); -[data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00'); -data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); -options.INTERP = 1; -options.POLARPLOT = 0; -options.BWR = 0; % bluewhitered plot or gray -options.CLIMAUTO = 1; % adjust the colormap auto -% options.NAME = '\phi'; -% options.NAME = '\phi^{NZ}'; -% options.NAME = '\omega_z'; -options.NAME = 'N_i^{00}'; -% options.NAME = 's_{Ey}'; -% options.NAME = 'n_i^{NZ}'; -% options.NAME = 'Q_x'; -% options.NAME = 'n_i'; -% options.NAME = 'n_i-n_e'; -options.PLAN = 'kxky'; -% options.NAME = 'f_i'; -% options.PLAN = 'sx'; -options.COMP = 'avg'; -% options.TIME = data.Ts5D(end-30:end); -options.TIME = data.Ts3D(1:1:end); -% options.TIME = [0:1500]; -data.EPS = 0.1; -data.a = data.EPS * 2000; -options.RESOLUTION = 256; -create_film(data,options,'.gif') -end - -if 0 -%% field snapshots +if 1 +%% 2D field snapshots % Options -[data.Na00, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'Na00'); -[data.PHI, data.Ts3D] = compile_results_3D(data.folder,J0,J1,'phi'); -data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); - options.INTERP = 0; options.POLARPLOT = 0; options.AXISEQUAL = 0; options.NORMALIZE = 0; options.LOGSCALE = 0; options.CLIMAUTO = 1; -options.NAME = 'N_i^{00}'; -% options.NAME = 's_{Ey}'; +options.TAVG = 1; +% options.NAME = ['N_e^{00}']; +% options.NAME = 'n_e'; +% % options.NAME = 'u_i'; +% options.NAME = 'T_i'; +% options.NAME = 'Q_{xe}'; +options.NAME = 'v_{Ey}'; +% options.NAME = 'w_{Ez}'; +% options.NAME = '\omega_z'; % options.NAME = '\phi'; -options.PLAN = 'yz'; -options.COMP = 'avg'; -options.TIME = [50 100]; -% options.TIME = data.Ts3D(1:2:end); +% options.NAME = 'n_i-n_e'; +loc =11; +[~,i_] = min(abs(loc - data.grids.y)); +options.COMP =i_; +% options.PLAN = '3D'; +options.PLAN = 'xy'; options.COMP =floor(data.grids.Nz/2)+1; +% options.PLAN = 'yz'; options.COMP ='avg'; +% options.COMP ='avg'; +options.XYZ =[-11 20 0]; +options.TIME = [90:150]; options.RESOLUTION = 256; fig = photomaton(data,options); -colormap(gray) +% colormap(gray) clim('auto') % save_figure(data,fig) end @@ -107,42 +187,147 @@ if 0 profiler(data) end -if 0 +if 1 +%% Mode evolution +% [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi'); +% [data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00'); +% data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +% data.Ne00 = reshape(data.Na00(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); + +options.NORMALIZED = 0; +options.TIME = data.Ts3D; +options.KX_TW = [0.1 2.5]; %kx Growth rate time window +options.KY_TW = [0.1 2.5]; %ky Growth rate time window +options.NMA = 1; +options.NMODES = 64; +options.iz = 'avg'; % avg or index +options.ik = 1; % sum, max or index +options.fftz.flag = 0; +options.FIELD = 'Ni00'; +% options.FIELD = 'phi'; +% options.FIELD = 'T_i'; +options.GOK2 = 0; +options.SHOWFIG = 1; +[fig, wkykx, ekykx] = mode_growth_meter(data,options); +%% +kx = (1:data.grids.Nx/2)'*2*pi/data.fort_00.GRID.Lx; +ky = (1:data.grids.Ny/2)'*2*pi/data.fort_00.GRID.Ly; +gkxky = real(wkykx(2:end,1:data.grids.Nx/2))'; +gkxky(isnan(gkxky)) =0; +gkxky(isinf(gkxky)) =0; +% gkxky(gkxky<0) =0; +% gkxky = imgaussfilt(gkxky,1); +% +wkxky = imag(wkykx(2:end,1:data.grids.Nx/2))'; +wkxky(isnan(wkxky)) =0; +wkxky(isinf(wkxky)) =0; +% wkxky(wkxky<0) =0; +% wkxky = imgaussfilt(wkxky,1.5); +% +figure; +subplot(121) + contourf(kx,ky,gkxky',10) + % clim(0.5*[0 1]); + % colormap(bluewhitered); colorbar; + xlim([0.025 1]); + xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$') +subplot(122) + contourf(kx,ky,wkxky',10) + % clim(1*[0 1]); + % colormap(bluewhitered); colorbar + xlim([0.025 1]); + xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$') +% save_figure(data,fig,'.png') +end + +if 1 %% Hermite-Laguerre spectrum -[data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Napjz'); +[data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,0,10,'Napjz'); % [data.Napjz, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Nipjz'); -options.ST = 0; +options.ST = 1; options.NORMALIZED = 0; -options.LOGSCALE = 1; +options.LOGSCALE = 0; options.FILTER = 0; %filter the 50% time-average of the spectrum from options.TAVG_2D = 0; %Show a 2D plot of the modes, 50% time averaged options.TAVG_2D_CTR= 0; %make it contour plot +options.TWINDOW = [6 20]; fig = show_moments_spectrum(data,options); end -if 0 -%% Mode evolution -[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi'); -[data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00'); -data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +if (0 && NZ>4) +%% Ballooning plot +% [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi'); +if data.inputs.BETA > 0 +[data.PSI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'psi'); +end +options.time_2_plot = [25 100]; +options.kymodes = 1.5; +options.normalized = 1; +options.PLOT_KP = 0; +% options.field = 'phi'; +options.SHOWFIG = 1; +[fig, chi, phib, psib, ~] = plot_ballooning(data,options); +end -options.NORMALIZED = 0; -options.TIME = data.Ts3D; -options.KX_TW = [ 0 20]; %kx Growth rate time window -options.KY_TW = [ 0 50]; %ky Growth rate time window -options.NMA = 1; -options.NMODES = 50; -options.iz = 'avg'; % avg or index -options.ik = 9; % sum, max or index -options.fftz.flag = 0; -% options.FIELD = 'Ni00'; -options.FIELD = 'phi'; -options.GOK2 = 0; -fig = mode_growth_meter(data,options); -% save_figure(data,fig,'.png') +if 1 +%% 1D spectral plot +options.TIME = [30 80]; % averaging time window +% options.NAME = ['N_i^{00}']; +% options.NAME = 'n_i'; +% options.NAME = 'T_i'; +% options.NAME = 'Q_{xi}'; +% options.NAME = 's_{Ey}'; +options.NAME = '\psi'; +options.NORMALIZE = 0; +[fig] = plot_spectrum(data,options); end +if 0 +%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Options +% [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi'); +% [data.PSI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'psi'); +% [data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00'); +% data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +% data.Ne00 = reshape(data.Na00(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +% [data.DENS, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'dens'); +% data.DENS_I = reshape(data.DENS(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +% [data.TEMP, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'temp'); +% data.TEMP_I = reshape(data.TEMP(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D)); +options.INTERP = 0; +options.POLARPLOT = 0; +options.BWR = 1; % bluewhitered plot or gray +options.CLIMAUTO = 0; % adjust the colormap auto +% options.NAME = '\phi'; +% options.NAME = 'w_{Ez}'; +% options.NAME = '\psi'; +options.NAME = 'T_i'; +% options.NAME = '\phi^{NZ}'; +% options.NAME = ['N_e^{00}']; +% options.NAME = ['N_i^{00}']; +options.PLAN = 'xy'; +% options.PLAN = '3D'; +% options.XYZ =[-11 20 0]; +% options.COMP = 'avg'; +options.COMP = floor(data.grids.Nz/2+1); +options.TIME = data.Ts3D(1:1:end); +% options.TIME = [0:1500]; +data.EPS = 0.1; +data.a = data.EPS * 2000; +options.RESOLUTION = 256; +options.FPS = 12; +options.RMAXIS = 0; +create_film(data,options,'.gif'); +end + +if 0 +%% Metric infos +options.SHOW_FLUXSURF = 1; +options.SHOW_METRICS = 1; +[fig, geo_arrays] = plot_metric(data,options); +end + if 0 %% Study singular values [data.SV_ky_pj, data.Ts2D] = compile_results_2D(DATADIR,J0,J1,'sv_ky_pj'); diff --git a/wk/lin_run_script.m b/wk/lin_run_script.m index f872fc8767ddeeb77a2475406efafb4050956861..7e0eb7271fab42b84a51e3fea2c12c986b5eee4b 100644 --- a/wk/lin_run_script.m +++ b/wk/lin_run_script.m @@ -27,15 +27,17 @@ EXECNAME = 'gyacomo23_dp'; % double precision % run lin_DTT_HM_rho98 % run lin_DIIID_LM_rho90 % run lin_DIIID_LM_rho95 +% run lin_DIIID_LM_rho95_HEL % run lin_JET_rho97 % run lin_Entropy % run lin_ITG % run lin_KBM % run lin_RHT -rho = 0.95; TRIANG = 'NT'; READPROF = 0; -prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/']; +% run lin_Ivanov +rho = 0.95; TRIANG = 'NT'; READPROF = 1; +% prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/']; % prof_folder = ['parameters/profiles/DIIID_Oak_Nelson/',TRIANG,'/']; -% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/']; +prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/']; run lin_DIIID_data % run lin_STEP_EC_HD_psi71 % run lin_STEP_EC_HD_psi49 @@ -43,47 +45,31 @@ if 1 % Plot the profiles plot_params_vs_rho(geom,prof_folder,rho,0.5,1.1,Lref,mref,Bref,READPROF); end +% SIMID = ['rho_scan_DIIID_AUSTIN_2019/3x2x192x96x32/rho',num2str(rho)]; %% Change parameters -% EXBRATE = 0.0; % Background ExB shear flow -% K_Ti = 5.3; -% NU = 0.001; -CO = 'LD'; -GKCO = 1; -kymin= 0.3; LY = 2*pi/kymin; +% GEOMETRY = 's-alpha'; +DELTA =0.0; +K_Ni = 0; K_Ne = 0; +% DELTA = 0.0; +% DELTA = 0.2; +S_DELTA = DELTA/2; +LY = 2*pi/0.25; +TMAX = 20; NY = 2; -NX = 4; -NZ = 32; -PMAX = 2; -JMAX = PMAX/2; -DT = 20e-4; -EXBRATE = 0; -% EPS = 0.25; -% KAPPA = 1.0; -S_DELTA = min(2.0,S_DELTA); -% DELTA = -DELTA; -% PT parameters -% TAU = 5.45; -% K_Ne = 0; %vs 2.79 -% K_Ni = 0; -% K_Te = 9.6455;%vs 17.3 -% K_Ti = 3.3640;%vs 5.15 -% BETA = 7.9e-4; -% NU = 0.1; -MU_X = 0.0; MU_Y = 0.0; -SIGMA_E = 0.04; -TMAX = 1; -% DTSAVE0D = 200*DT; -DTSAVE3D = 0.1; -%%------------------------------------------------------------------------- +DT = 0.01; +TAU = 1; NU = 0.05; +% TAU = 1e-3; K_Ti = K_Ti/2/TAU; NU = 3*NU/8/TAU; ADIAB_E = 1; NA = 1; +% MU_X = 1; MU_Y = 1; %% RUN setup +system(['cp ../results/',SIMID,'/',PARAMS,'/fort_00.90 ../results/',SIMID,'/.']); % system(['rm fort*.90']); % Run linear simulation if RUN MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;']; % RUN =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;']; - RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;']; - % RUN =['time ',mpirun,' -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0;']; + RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0;']; + % RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 2 2 1 0;']; % RUN =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;']; % RUN =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;']; % RUN = ['./../../../bin/gyacomo23_sp 0;']; @@ -102,15 +88,15 @@ data = {}; % Initialize data as an empty cell array % load grids, inputs, and time traces data = compile_results_low_mem(data,LOCALDIR,J0,J1); -if 0 % Activate or not +if 1 % Activate or not %% plot mode evolution and growth rates % Load phi [data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi'); options.NORMALIZED = 0; options.TIME = data.Ts3D; % Time window to measure the growth of kx/ky modes -options.KY_TW = [0.5 1.0]*data.Ts3D(end); -options.KX_TW = [0.5 1.0]*data.Ts3D(end); +options.KY_TW = [0.25 1.0]*data.Ts3D(end); +options.KX_TW = [0.25 1.0]*data.Ts3D(end); options.NMA = 1; % Set NMA option to 1 options.NMODES = 999; % Set how much modes we study options.iz = 'avg'; % Compressing z @@ -118,8 +104,38 @@ options.ik = 1; % options.GOK2 = 0; % plot gamma/k^2 options.fftz.flag = 0; % Set fftz.flag option to 0 options.FIELD = 'phi'; -options.SHOWFIG = 1; -[fig, kykx, wkykx, ekykx] = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig +options.SHOWFIG = 1; +[fig, wkykx, ekykx] = mode_growth_meter(data,options); +% %% +% kx = (1:data.grids.Nx/2)'*2*pi/data.fort_00.GRID.Lx; +ky = (1:data.grids.Ny/2)'*2*pi/data.fort_00.GRID.Ly; +gkxky = real(wkykx(2:end,1:data.grids.Nx/2))'; +gkxky(isnan(gkxky)) =0; +gkxky(isinf(gkxky)) =0; +figure; plot(ky,gkxky(1,:)); +% gkxky(gkxky<0) =0; +% % gkxky = imgaussfilt(gkxky,1); +% % +% wkxky = imag(wkykx(2:end,1:data.grids.Nx/2))'; +% wkxky(isnan(wkxky)) =0; +% wkxky(isinf(wkxky)) =0; +% % wkxky(wkxky<0) =0; +% % wkxky = imgaussfilt(wkxky,1.5); +% % +% figure; +% subplot(121) +% contourf(kx,ky,gkxky',10) +% % clim(0.5*[0 1]); +% % colormap(bluewhitered); colorbar; +% xlim([0.025 1]); +% xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$') +% subplot(122) +% contourf(kx,ky,wkxky',10) +% % clim(1*[0 1]); +% % colormap(bluewhitered); colorbar +% xlim([0.025 1]); +% xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$') +% % save_figure(data,fig,'.png') end if (0 && NZ>4) @@ -144,13 +160,13 @@ ikx = 2; iky = 1; t0 = 1; 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)); t_ = data.Ts3D(it0:it1); -% TH = 1.635*EPS^1.5 + 0.5*EPS^2+0.35*EPS^2.5; theory = 1/(1+Q0^2*TH/EPS^2); +TH = 1.635*EPS^1.5 + 0.5*EPS^2+0.36*EPS^2.5; theory = 1/(1+Q0^2*TH/EPS^2); clr_ = lines(20); figure -plot(t_, plt(data.PHI)); hold on; +plot(t_, -plt(data.PHI)); hold on; plot(t_,0.5* exp(-t_*NU)+theory,'--k'); plot([t_(1) t_(end)],theory*[1 1],'-k'); -plot([t_(1) t_(end)],mean(plt(data.PHI))*[1 1],'-g'); +plot([t_(1) t_(end)],-mean(plt(data.PHI))*[1 1],'-g'); xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$') title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.grids.kx(ikx),data.grids.ky(iky))) end @@ -160,7 +176,7 @@ if 0 plot_metric(data); end -if 1 +if 0 %% Compiled plot for lin EM analysis [data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi'); if data.inputs.BETA > 0 @@ -185,16 +201,22 @@ if 1 options.FIELD = 'phi'; options.SHOWFIG = 0; [~, chi, phib, psib, ~] = plot_ballooning(data,options); - [~, kykx, wkykx, ekykx] = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig + [fig, wkykx, ekykx] = mode_growth_meter(data,options); + kx = (1:data.grids.Nx/2)'*2*pi/data.fort_00.GRID.Lx; + ky = (1:data.grids.Ny/2)'*2*pi/data.fort_00.GRID.Ly; [~,~,R,Z] = plot_metric(data,options); figure; subplot(3,2,1); plot(chi,real(phib),'-b'); hold on; plot(chi,imag(phib),'-r'); xlabel('$\chi$'); ylabel('$\phi$') subplot(3,2,3); plot(chi,real(psib),'-b'); hold on; plot(chi,imag(psib),'-r'); xlabel('$\chi$'); ylabel('$\psi$') - subplot(3,2,5); plot(squeeze(kykx(:,1)),squeeze(real(wkykx(:,1))),'o--'); hold on; - plot(squeeze(kykx(:,1)),squeeze(imag(wkykx(:,1))),'o--'); - xlabel('$k_y\rho_s$'); ylabel('$\gamma,\omega$'); + subplot(3,2,5); errorbar(ky,squeeze(real(wkykx(2:end,1))),... + squeeze(real(ekykx(2:end,1))),... + 'ok--','MarkerSize',8,'LineWidth',1.5); hold on; + errorbar(ky,squeeze(imag(wkykx(2:end,1))),... + squeeze(imag(ekykx(2:end,1))),... + '*k--','MarkerSize',8,'LineWidth',1.5); + xlabel('$k_y\rho_s$'); ylabel('$\gamma (o),\,\omega (*)$'); R = R*geom.R0; Z = Z*geom.R0; subplot(1,2,2); plot(R,Z,'-k'); xlabel('$R$'); ylabel('$Z$'); axis equal diff --git a/wk/lin_scan_script.m b/wk/lin_scan_script.m index bf36d44d2d156a2c46e53834f47edaa055c7b349..3883f0e9b4fff3be942ca2dc74e8f14c0b89dbbc 100644 --- a/wk/lin_scan_script.m +++ b/wk/lin_scan_script.m @@ -28,24 +28,29 @@ EXECNAME = 'gyacomo23_dp'; % double precision % run lin_Entropy % run lin_ITG % run lin_RHT -rho = 0.95; TRIANG = 'NT'; READPROF = 1; +rho = 0.95; TRIANG = 'PT'; READPROF = 1; % prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/']; % prof_folder = ['parameters/profiles/DIIID_Oak_Nelson/',TRIANG,'/']; prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/']; run lin_DIIID_data %% Change parameters -NU = 1; -TAU = 1; +% NU = 1; +% TAU = 1; NY = 2; EXBRATE = 0; +S_DELTA = min(2.0,S_DELTA); SIGMA_E = 0.023; +NEXC = 0; +LX = 120; %% Scan parameters SIMID = [SIMID,'_scan']; P_a = [2 4]; % P_a = 2; -ky_a = [0.01 0.02 0.05 0.1 0.2 0.5 1.0 2.0 5.0 10.0]; -CO = 'LD'; +ky_a = [0.01 0.02 0.05 0.1 0.2 0.5 1.0 2.0 5.0 10.0]; +% ky_a = 4.0; +dt_a = logspace(-2,-3,numel(ky_a)); +CO = 'DG'; %% Scan loop % arrays for the result g_ky = zeros(numel(ky_a),numel(P_a)); @@ -58,10 +63,10 @@ for PMAX = P_a i = 1; for ky = ky_a LY = 2*pi/ky; - DT = 2e-4;%/(1+log(ky/0.05));%min(1e-2,1e-3/ky); - TMAX = 20;%min(10,1.5/ky); - DTSAVE0D = 0.1; - DTSAVE3D = 0.01; + DT = dt_a(i);%1e-3;%/(1+log(ky/0.05));%min(1e-2,1e-3/ky); + TMAX = DT*10000;%2;%min(10,1.5/ky); + DTSAVE0D = 100*DT; + DTSAVE3D = 10*DT; %% RUN setup % naming @@ -94,21 +99,27 @@ for PMAX = P_a data_ = compile_results_low_mem(data_,LOCALDIR,00,00); [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi'); - - % linear growth rate (adapted for 2D zpinch and fluxtube) - options.TRANGE = [0.5 1]*data_.Ts3D(end); - options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z - options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3 - - [~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window - [~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ... - [wkykx,ekykx] = compute_growth_rates(data_.PHI(:,:,:,it1:it2),data_.Ts3D(it1:it2)); + options.NORMALIZED = 0; + options.TIME = data_.Ts3D; + % Time window to measure the growth of kx/ky modes + options.KY_TW = [0.7 1.0]*data_.Ts3D(end); + options.KX_TW = [0.7 1.0]*data_.Ts3D(end); + options.NMA = 1; % Set NMA option to 1 + options.NMODES = 999; % Set how much modes we study + options.iz = 'avg'; % Compressing z + options.ik = 1; % + options.GOK2 = 0; % plot gamma/k^2 + options.fftz.flag = 0; % Set fftz.flag option to 0 + options.FIELD = 'phi'; + options.SHOWFIG = 0; + [fig, wkykx, ekykx] = mode_growth_meter(data_,options); + % [wkykx,ekykx] = compute_growth_rates(data_.PHI(:,:,:,it1:it2),data_.Ts3D(it1:it2)); g_ky (i,j) = real(wkykx(2,1)); g_std(i,j) = real(ekykx(2,1)); w_ky (i,j) = imag(wkykx(2,1)); w_std(i,j) = imag(ekykx(2,1)); [gmax, ikmax] = max(g_ky(i,j)); - + msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg); end i = i + 1; diff --git a/wk/load_metadata_scan.m b/wk/load_metadata_scan.m index a9df5c95be4a524979dbae0448bf4f1f2461c3c2..89db4acbe3035dde9e6d5f3226eba85105138ded 100644 --- a/wk/load_metadata_scan.m +++ b/wk/load_metadata_scan.m @@ -11,7 +11,7 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0' % datafname = 'lin_DIIID_Oak_Nelson_high_density_NT_scan/6x32_ky_0.01_10_P_2_8_kN_1.0883_LDGK_0.0080915_be_0.0015991.mat'; % rho = 0.95 % datagname = 'lin_DIIID_Oak_Nelson_high_density_PT_scan/6x32_ky_0.01_10_P_2_8_kN_0.62888_LDGK_0.0046858_be_0.0048708.mat'; -datafname = 'lin_DIIID_Oak_Nelson_high_density_NT_scan/6x32_ky_0.01_10_P_2_8_kN_1.6989_LDGK_0.0167_be_0.00075879.mat'; +datafname = 'lin_DIIID_Oak_Nelson_high_density_PT_scan/6x32_ky_0.01_10_P_2_4_kN_0.62888_DGGK_0.0046858_be_0.0048708.mat'; %% Chose if we filter gamma>0.05 FILTERGAMMA = 1; diff --git a/wk/old scripts/p3_geom_scan_analysis.m b/wk/old scripts/p3_geom_scan_analysis.m new file mode 100644 index 0000000000000000000000000000000000000000..1c63315622f6a2a84a842307f9f66e86d04c9a7f --- /dev/null +++ b/wk/old scripts/p3_geom_scan_analysis.m @@ -0,0 +1,306 @@ +% Get the current directory (wk) +curdir = pwd; +NCONTOUR = 0; +% give ref value and param names +REFVAL= 0; +% normalize to the max all data +NORM_ALL= 0; +% normalize to the max each line +NORM_LIN= 0; +% normalize to the max each column +NORM_COL= 0; +% Get and plot the fluxsurface +GETFLUXSURFACE = 0; + +% partition= '../results/paper_3/'; +% Get the scan directory +switch 6 + case 1 % delta kappa scan + casename = 'DTT rho85'; + partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/'; + scandir = 'P2_J1_delta_kappa_scan'; scanname= '(2,1)'; + % scandir = 'P4_J2_delta_kappa_scan'; scanname= '(4,2)'; + nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23; scale1 =1.0; + nml2 = 'GEOMETRY'; pnam2 = '$\kappa$'; attr2 = 'kappa'; pref2 = 1.53; scale2 =1.0; + t1 = 50; t2 = 150; + case 2 % shear safety factor scan + casename = 'DTT rho85'; + partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/'; + scandir = 'P2_J1_PT_sfact_shear_scan'; scanname= '(2,1)'; + % scandir = 'P2_J1_NT_sfact_shear_scan'; scanname= '(2,1)'; + nml1 = 'GEOMETRY'; pnam1 = '$\hat s$'; attr1 = 'shear'; pref1 = 3.63; scale1 =1.0; + nml2 = 'GEOMETRY'; pnam2 = '$q_0$'; attr2 = 'q0'; pref2 =-2.15; scale2 =1.0; + t1 = 50; t2 = 150; + case 3 + casename = 'DTT rho85'; + partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/'; + % scandir = 'P2_J1_delta_nuDGGK_scan'; scanname= 'DG (2,1)'; + % scandir = 'P4_J2_delta_nuDGGK_scan'; scanname= 'DG (4,2)'; + % scandir = 'P8_J4_delta_nuDGGK_conv_test'; scanname= 'DG (8,4)'; + % scandir = 'P2_J1_delta_nuSGGK_scan'; scanname= 'SG (2,1)'; + % scandir = 'P4_J2_delta_nuSGGK_scan'; scanname= 'SG (4,2)'; + % scandir = 'P8_J4_delta_nuSGGK_conv_test'; scanname= 'SG (8,4)'; + % scandir = 'P4_J2_delta_nuSGGKii_scan'; scanname= 'SGii (4,2)'; + scandir = 'P2_J1_delta_nuLDGK_scan'; scanname= 'LD (2,1)'; + % scandir = 'P4_J2_delta_nuLDGK_scan'; scanname= 'LD (4,2)'; + nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23; scale1 =1.0; + nml2 = 'MODEL'; pnam2 = '$\nu$'; attr2 = 'nu'; pref2 = 0.5; scale2 =1.0; + t1 = 50; t2 = 150; + case 4 % shear delta scan + casename = 'DIIID rho95'; + partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/'; + scandir = '3x2x192x48x32_RT_2500_delta_shear_scan'; scanname= 'CI DG RT 2500'; + % scandir = '3x2x256x64x48_delta_shear_scan'; + nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0; + nml2 = 'GEOMETRY'; pnam2 = '$\hat s$'; attr2 = 'shear'; pref2 = 0.8; scale2 =1.0; + t1 = 50; t2 = 150; + case 5 % delta K_T tau=1 + casename = 'DIIID rho95 $\tau=1$'; + partition= '/misc/gyacomo23_outputs/paper_3/DIIID_tau_1_rho95_geom_scan/'; + scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)'; + % scandir = '5x3x192x48x32_delta_RT_scan'; scanname= '(4,2)'; + % scandir = 'delta_RT_scan_PJ_21'; scanname= '(2,1)'; + nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0; + nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1.0; + t1 = 200; t2 = 500; + case 6 % delta K_T cold ions + casename = 'DIIID rho95 $\tau=10^{-3}$'; + partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/'; + scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)'; + nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0; + nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1e3/2; + t1 = 200; t2 = 480; + case 7 % delta s_delta + casename = 'DIIID rho95 $\tau=10^{-3}$'; + partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/'; + scandir = '3x2x192x48x32_RT_1000_delta_sdelta_scan'; scanname= 'RT=1000 (2,1)'; + % scandir = ''; + nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0; + nml2 = 'GEOMETRY'; pnam2 = '$s_\delta$'; attr2 = 's_delta'; pref2 = 0.8; scale2 =1.0; + t1 = 200; t2 = 295; + case 8 % eps q0 + casename = 'DIIID rho95 $\tau=10^{-3}$'; + partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/'; + scandir = '3x2x192x48x32_RT_1000_eps_q0_scan/PT'; scanname= 'PT, RT=1000 (2,1)'; + % scandir = '3x2x192x48x32_RT_1000_eps_q0_scan/NT'; scanname= 'NT, RT=1000 (2,1)'; + % scandir = ''; + nml1 = 'GEOMETRY'; pnam1 = '$\epsilon$'; attr1 = 'eps'; pref1 = 0; scale1 =1.0; + nml2 = 'GEOMETRY'; pnam2 = '$q_0$'; attr2 = 'q0'; pref2 = 0; scale2 =1.0; + t1 = 200; t2 = 400; + case 9 % CBC Dimits shift + casename = 'HEL CBC'; + partition= '/misc/gyacomo23_outputs/paper_3/HEL_CBC/'; + scandir = '128x32x24'; scanname= 'CBC HEL'; + nml1 = 'SPECIES'; pnam1 = '$R_T$'; attr1 = 'k_T_'; pref1 = 0; scale1 =1.0; + nml2 = 'SPECIES'; pnam2 = '$R_N$'; attr2 = 'k_N_'; pref2 = 0; scale2 =1.0; + t1 = 1000; t2 = 2000; +end +scanname= [casename scanname]; +scandir = [partition,scandir,'/']; +% Get a list of all items in the current directory +contents = dir(scandir); + +% Iterate through the contents +Qxavg = []; Qxerr = []; para1 = []; para2 = []; R = []; Z = []; +Qxt = struct(); +for i = 1:length(contents) + % Check if the item is a directory and not '.' or '..' + if contents(i).isdir && ~strcmp(contents(i).name, '.') ... + && ~strcmp(contents(i).name, '..') + % Get and display the name of the subdirectory + subdir = [scandir,contents(i).name]; + disp(['Subdirectory: ' contents(i).name]); + % Get parameters + param = read_namelist([subdir,'/fort_00.90']); + para1 = [para1 param.(nml1).(attr1)]; + para2 = [para2 param.(nml2).(attr2)]; + % Now you are in the subdirectory. You can perform operations here. + [t_all, Pxi_all, Qxi_all, Pxe_all, Qxe_all] = read_flux_out_XX(subdir); + if(numel(Qxe_all) > 1) + Qxtot = Qxi_all+Qxe_all; + else + Qxtot = Qxi_all; + end + Qxt.(['dat_',num2str(i)]) = struct(); + Qxt.(['dat_',num2str(i)]).Qx = Qxtot; + Qxt.(['dat_',num2str(i)]).t = t_all; + Qxt.(['dat_',num2str(i)]).name = contents(i).name; + if(numel(t_all) > 1) + disp(num2str(t_all(end))) + [~,it1] = min(abs(t_all-t1)); + [~,it2] = min(abs(t_all-t2)); + steady_slice = it1:it2; + if(t_all(end) >= t2) + [fullAvg,sliceAvg,sliceErr] = sliceAverage(Qxtot(steady_slice),3); + Qxavg = [Qxavg fullAvg]; + Qxerr = [Qxerr mean(sliceErr)]; + else + Qxavg = [Qxavg nan]; + Qxerr = [Qxerr nan]; + end + else + Qxavg = [Qxavg nan]; + Qxerr = [Qxerr nan]; + end + end + if GETFLUXSURFACE + data = load([subdir,'/RZ.txt']); + R_ = data(:, 1); + Z_ = data(:, 2); + R_ = [R_;R_(1)]'; Z_ = [Z_;Z_(1)]'; + R = [R ; R_]; Z = [Z ; Z_]; + end + +end +if 0 +%% plot time traces +attr = fieldnames(Qxt); +Nsim = numel(attr); +figure +% compute growth at the begining +tw = [10 40]; +gr = 1:Nsim; err = 1:Nsim; +for i = 1:1:Nsim + tmp_ = Qxt.(attr{i}); + t = tmp_.t; + y = tmp_.Qx; + plot(t,y,'DisplayName',tmp_.name); hold on; + [~,it1] = min(abs(t-tw(1))); + [~,it2] = min(abs(t-tw(2))); + [gr_, err_] = compute_growth(t(it1:it2),y(it1:it2)); + gr(i) = gr_; err(i) = err_; +end +%% +toplot = real(reshape(gr,sz))'; +toplot = toplot(idx1,idx2); + +figure +imagesc_custom(xx_,yy_,toplot); hold on +end +%% reshaping, sorting and plotting +p1 = unique(para1)/scale1; +p2 = unique(para2)/scale2; +N1 = numel(p1); +N2 = numel(p2); + +if para1(1) == para1(2) + sz = [N2 N1]; + TRANSPOSE = 1; +else + sz = [N1 N2]; + TRANSPOSE = 0; +end + +Zavg = reshape(Qxavg,sz); +Zerr = reshape(Qxerr,sz); +XX = reshape(para1/scale1,sz); +YY = reshape(para2/scale2,sz); + +if TRANSPOSE + Zavg = Zavg'; + Zerr = Zerr'; + XX = XX'; + YY = YY'; +end + +[~,idx1] = sort(XX(:,1)); +[~,idx2] = sort(YY(1,:)); +Zavg = Zavg(idx1,idx2); +Zerr = Zerr(idx1,idx2); +XX = XX(idx1,idx2); +YY = YY(idx1,idx2); + +% compute the +if REFVAL + if NORM_ALL + Qxname = '$\bar Q_{tot}/\bar Q_{max}[\%]$'; + [tmp,iref1] = max(Zavg); + [~, iref2] = max(tmp); + iref1 = iref1(iref2); + else + Qxname = '$\langle (Q_{tot}-Q_{ref})/Q_{ref} \rangle_t[\%]$'; + if pref1 ~= 999 + [~,iref1] = min(abs(XX(:,1)-pref1)); + else + iref1 = 1:N1; + end + if pref2 ~= 999 + [~,iref2] = min(abs(YY(1,:)-pref2)); + else + iref2 = 1:N2; + end + end + iref1 = ones(N1,1).*iref1; + iref2 = ones(N2,1).*iref2; + xref = XX(iref1,iref2); + yref = YY(iref1,iref2); + Qxref = Zavg(iref1,iref2); + Qrefname = ['$Q_{ref}=$',num2str(Qxref(1,1))]; +else + Qref = 1; + if NORM_LIN + Qxname = '$\bar Q_{tot}/\bar Q_{max}[\%]$, per line'; + for il = 1:sz(1) + maxline = max(Zavg(:,il)); + Zavg(:,il) = Zavg(:,il)./maxline; + Zerr(:,il) = Zerr(:,il)./maxline; + end + elseif NORM_COL + Qxname = '$\bar Q_{tot}/\bar Q_{max}[\%]$, per column'; + for ic = 1:sz(2) + maxcol = max(Zavg(ic,:)); + Zavg(ic,:) = Zavg(ic,:)./maxcol; + Zerr(ic,:) = Zerr(ic,:)./maxcol; + end + else + Qxname = '$\langle Q_{tot} \rangle_t$'; + end +end + +% Figure +figure +subplot(1,2,1) +[xx_,yy_] = meshgrid(XX(:,1),YY(1,:)); +if REFVAL + if NORM_ALL || NORM_COL || NORM_LIN + toplot = (Zavg./Qxref)' + CLIM = [0 1]; + else + toplot = ((Zavg-Qxref)./Qxref * 100)'; + CLIM = 'auto'; + end +else + toplot = Zavg'; + CLIM = 'auto'; +end +if NCONTOUR <= 0 + imagesc_custom(xx_,yy_,toplot); hold on +else + contourf(XX(:,1),YY(1,:),Zavg',NCONTOUR); hold on +end +if REFVAL && ~((pref1==999) || (pref2==999)) + plot(xref(1,1),yref(1,1),'xk','MarkerSize',14,'DisplayName',Qrefname) + legend('show') +end +xlabel(pnam1); ylabel(pnam2); +title(scanname) +colormap(bluewhitered); colorbar; clim(CLIM); +if ~REFVAL + colormap(jet); +end +subplot(1,2,2) +clrs = jet(N2); +for i = 1:N2 + errorbar(XX(:,i),Zavg(:,i),Zerr(:,i),... + 'DisplayName',[pnam2,'=',num2str(p2(i))],... + 'Color',clrs(i,:)); + hold on; +end +if REFVAL && ~((pref1==999) || (pref2==999)) + plot(xref(1,1),0,'xk','MarkerSize',14,'DisplayName',Qrefname) +end +grid on +xlabel(pnam1); ylabel('$\langle Q_{tot} \rangle_t$'); +legend('show','Location','northwest'); +title([param.COLLISION.collision_model{1}, ... + ', $(P,J)=(',num2str(param.GRID.pmax),',',num2str(param.GRID.jmax),')$']) diff --git a/wk/p3_geom_scan_analysis.m b/wk/p3_geom_scan_analysis.m index 0ae2554962bbe6fe28fa0ce57bb36140f3c165a3..d81207835636db7190a19f13bcbfc7d9d95f7581 100644 --- a/wk/p3_geom_scan_analysis.m +++ b/wk/p3_geom_scan_analysis.m @@ -14,87 +14,38 @@ GETFLUXSURFACE = 0; % partition= '../results/paper_3/'; % Get the scan directory -switch 6 - case 1 % delta kappa scan - casename = 'DTT rho85'; - partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/'; - scandir = 'P2_J1_delta_kappa_scan'; scanname= '(2,1)'; - % scandir = 'P4_J2_delta_kappa_scan'; scanname= '(4,2)'; - nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23; scale1 =1.0; - nml2 = 'GEOMETRY'; pnam2 = '$\kappa$'; attr2 = 'kappa'; pref2 = 1.53; scale2 =1.0; - t1 = 50; t2 = 150; - case 2 % shear safety factor scan - casename = 'DTT rho85'; - partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/'; - scandir = 'P2_J1_PT_sfact_shear_scan'; scanname= '(2,1)'; - % scandir = 'P2_J1_NT_sfact_shear_scan'; scanname= '(2,1)'; - nml1 = 'GEOMETRY'; pnam1 = '$\hat s$'; attr1 = 'shear'; pref1 = 3.63; scale1 =1.0; - nml2 = 'GEOMETRY'; pnam2 = '$q_0$'; attr2 = 'q0'; pref2 =-2.15; scale2 =1.0; - t1 = 50; t2 = 150; - case 3 - casename = 'DTT rho85'; - partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/'; - % scandir = 'P2_J1_delta_nuDGGK_scan'; scanname= 'DG (2,1)'; - % scandir = 'P4_J2_delta_nuDGGK_scan'; scanname= 'DG (4,2)'; - % scandir = 'P8_J4_delta_nuDGGK_conv_test'; scanname= 'DG (8,4)'; - % scandir = 'P2_J1_delta_nuSGGK_scan'; scanname= 'SG (2,1)'; - % scandir = 'P4_J2_delta_nuSGGK_scan'; scanname= 'SG (4,2)'; - % scandir = 'P8_J4_delta_nuSGGK_conv_test'; scanname= 'SG (8,4)'; - % scandir = 'P4_J2_delta_nuSGGKii_scan'; scanname= 'SGii (4,2)'; - scandir = 'P2_J1_delta_nuLDGK_scan'; scanname= 'LD (2,1)'; - % scandir = 'P4_J2_delta_nuLDGK_scan'; scanname= 'LD (4,2)'; - nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23; scale1 =1.0; - nml2 = 'MODEL'; pnam2 = '$\nu$'; attr2 = 'nu'; pref2 = 0.5; scale2 =1.0; - t1 = 50; t2 = 150; - case 4 % shear delta scan - casename = 'DIIID rho95'; - partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/'; - scandir = '3x2x192x48x32_RT_2500_delta_shear_scan'; scanname= 'CI DG RT 2500'; - % scandir = '3x2x256x64x48_delta_shear_scan'; - nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0; - nml2 = 'GEOMETRY'; pnam2 = '$\hat s$'; attr2 = 'shear'; pref2 = 0.8; scale2 =1.0; - t1 = 50; t2 = 150; - case 5 % delta K_T tau=1 +switch 2 + case 1 % delta K_T tau=1 casename = 'DIIID rho95 $\tau=1$'; partition= '../results/paper_3/DIIID_tau_1_rho95_geom_scan/'; - % scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)'; + % % scandir = '3x2x192x48x32_nu_0.05_delta_RT_scan'; scanname= '(2,1)'; + % scandir = '3x2x192x48x32_nu_0.1_delta_RT_scan'; scanname= '(2,1)'; + % scandir = '3x2x192x48x24_nu_0.1_delta_RT_scan'; scanname= '(2,1)'; + % scandir = '3x2x192x48x32_nu_1.0_delta_RT_scan'; scanname= '(2,1)'; scandir = '2_1_delta_RT_scan'; scanname= '(2,1)'; - % scandir = '5x3x192x48x32_delta_RT_scan'; scanname= '(4,2)'; + % scandir = '5x3x192x48x32_nu_0.05_delta_RT_scan'; scanname= '(4,2)'; + scandir = '5x3x192x48x32_nu_1.0_delta_RT_scan'; scanname= '(4,2)'; + % scandir = '5x3x192x48x32_delta_RT_scan'; scanname= '(2,1)'; % scandir = 'delta_RT_scan_PJ_21'; scanname= '(2,1)'; nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0; - nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1.0; - t1 = 200; t2 = 500; - case 6 % delta K_T cold ions + nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 5; scale2 =1.0; + t1 = 300; t2 = 500; zfactor = 1; + case 2 % delta K_T cold ions casename = 'DIIID rho95 $\tau=10^{-3}$'; - partition= '../results/paper_3/DIIID_cold_ions_rho95_geom_scan/'; - scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)'; - nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0; - nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1e3/2; - t1 = 150; t2 = 280; - case 7 % delta s_delta - casename = 'DIIID rho95 $\tau=10^{-3}$'; - partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/'; - scandir = '3x2x192x48x32_RT_1000_delta_sdelta_scan'; scanname= 'RT=1000 (2,1)'; - % scandir = ''; + partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/'; + % scandir = '3x2x192x48x32_nu_0_delta_RT_scan'; scanname= '(2,1)'; + scandir = '3x2x192x48x32_nu_0.05_delta_RT_scan'; scanname= '(2,1)'; + nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1; + nml2 = 'SPECIES'; pnam2 = '$\kappa_T$'; attr2 = 'K_T_'; pref2 = 5; scale2 =500; + t1 = 80; t2 = 400; zfactor = 2; + case 3 % delta K_T HEL, better resolution + casename = 'DIIID rho95 $\tau=1$'; + % partition= '/misc/gyacomo23_outputs/paper_3/geom_scan_DIIID_HEL/NU_50/'; + partition= '/misc/gyacomo23_outputs/paper_3/geom_scan_DIIID_HEL/NU_20/'; + scandir = '.'; scanname= 'CBC HEL'; nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0; - nml2 = 'GEOMETRY'; pnam2 = '$s_\delta$'; attr2 = 's_delta'; pref2 = 0.8; scale2 =1.0; - t1 = 200; t2 = 295; - case 8 % eps q0 - casename = 'DIIID rho95 $\tau=10^{-3}$'; - partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/'; - scandir = '3x2x192x48x32_RT_1000_eps_q0_scan/PT'; scanname= 'PT, RT=1000 (2,1)'; - % scandir = '3x2x192x48x32_RT_1000_eps_q0_scan/NT'; scanname= 'NT, RT=1000 (2,1)'; - % scandir = ''; - nml1 = 'GEOMETRY'; pnam1 = '$\epsilon$'; attr1 = 'eps'; pref1 = 0; scale1 =1.0; - nml2 = 'GEOMETRY'; pnam2 = '$q_0$'; attr2 = 'q0'; pref2 = 0; scale2 =1.0; - t1 = 200; t2 = 400; - case 9 % CBC Dimits shift - casename = 'HEL CBC'; - partition= '/misc/gyacomo23_outputs/paper_3/HEL_CBC/'; - scandir = '128x32x24'; scanname= 'CBC HEL'; - nml1 = 'SPECIES'; pnam1 = '$R_T$'; attr1 = 'k_T_'; pref1 = 0; scale1 =1.0; - nml2 = 'SPECIES'; pnam2 = '$R_N$'; attr2 = 'k_N_'; pref2 = 0; scale2 =1.0; - t1 = 1000; t2 = 2000; + nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 5; scale2 =500; + t1 = 100; t2 = 150; zfactor = 1; end scanname= [casename scanname]; scandir = [partition,scandir,'/']; @@ -116,11 +67,16 @@ for i = 1:length(contents) para1 = [para1 param.(nml1).(attr1)]; para2 = [para2 param.(nml2).(attr2)]; % Now you are in the subdirectory. You can perform operations here. - [t_all, Pxi_all, Qxi_all, Pxe_all, Qxe_all] = read_flux_out_XX(subdir); + out = read_flux_out_XX(subdir); + t_all = out.t; + Pxi_all = out.Pxi; + Qxi_all = out.Qxi; + Pxe_all = out.Pxe; + Qxe_all = out.Qxe; if(numel(Qxe_all) > 1) - Qxtot = Qxi_all+Qxe_all; + Qxtot = zfactor*(Qxi_all+Qxe_all); else - Qxtot = Qxi_all; + Qxtot = zfactor*(Qxi_all); end Qxt.(['dat_',num2str(i)]) = struct(); Qxt.(['dat_',num2str(i)]).Qx = Qxtot; @@ -140,7 +96,7 @@ for i = 1:length(contents) Qxerr = [Qxerr nan]; end else - Qxavg = [Qxavg nan]; + Qxavg = [Qxavg nan]; Qxerr = [Qxerr nan]; end end @@ -159,7 +115,7 @@ attr = fieldnames(Qxt); Nsim = numel(attr); figure % compute growth at the begining -tw = [10 40]; +tw = [5 20]; gr = 1:Nsim; err = 1:Nsim; for i = 1:1:Nsim tmp_ = Qxt.(attr{i}); @@ -277,7 +233,7 @@ end if NCONTOUR <= 0 imagesc_custom(xx_,yy_,toplot); hold on else - contourf(XX(:,1),YY(1,:),Zavg',NCONTOUR); hold on + contour(XX(:,1),YY(1,:),Zavg'); hold on end if REFVAL && ~((pref1==999) || (pref2==999)) plot(xref(1,1),yref(1,1),'xk','MarkerSize',14,'DisplayName',Qrefname) @@ -305,3 +261,40 @@ xlabel(pnam1); ylabel('$\langle Q_{tot} \rangle_t$'); legend('show','Location','northwest'); title([param.COLLISION.collision_model{1}, ... ', $(P,J)=(',num2str(param.GRID.pmax),',',num2str(param.GRID.jmax),')$']) + +if 0 +%% plot minimum +idxmax = 1:numel(Zavg(1,:)); +idxmin = 1:numel(Zavg(1,:)); +xmax = 1:numel(Zavg(1,:)); +xmin = 1:numel(Zavg(1,:)); +ymax = 1:numel(Zavg(1,:)); +ymin = 1:numel(Zavg(1,:)); +err = 1:numel(Zavg(1,:)); + +x = linspace(min(p1),max(p1),128); +for i=1:numel(Zavg(1,:)) + [fit, dat] = polyfit(p1,Zavg(:,i)+0*Zerr(:,i),2); + [ymax(i),idx] = min(polyval(fit,x)); + xmax(i) = x(idx); + [fit, dat] = polyfit(p1,Zavg(:,i)-0*Zerr(:,i),2); + [ymin(i),idx] = min(polyval(fit,x)); + xmin(i) = x(idx); + + [zmin,idx] = min(Zavg(:,i)); + % err(i) = abs(zmin-Zavg(idx+1,i))/abs(zmin)+abs(zmin-Zavg(idx-1,i))/abs(zmin); +end +err = min(err,1); +xavg = 0.5*(xmax+xmin); +xerr = 0.5*abs(xmax-xmin); + +fit = polyfit(p2,xavg,1); +y = linspace(min(p2),max(p2),128); + +figure +plot(xavg+xerr,p2); hold on +plot(xavg-xerr,p2); hold on +plot(polyval(fit,y),y) +plot(polyval(fit,y),y) +plot(polyval(fit,y),y) +end \ No newline at end of file diff --git a/wk/parallel_scaling.m b/wk/parallel_scaling.m index ed35a92c497eedd0683c4d59d9a1d50616a051cb..5df3ebe177b84e9856bb87daadc5ea557ab760d8 100644 --- a/wk/parallel_scaling.m +++ b/wk/parallel_scaling.m @@ -1,11 +1,12 @@ -DIR = '/misc/gyacomo23_outputs/scaling/strong/17x9x256x256x32/'; scaling='strong'; +% DIR = '/misc/gyacomo23_outputs/scaling/strong/17x9x256x256x32/'; scaling='strong'; % DIR = '/misc/gyacomo23_outputs/scaling/strong/9x5x768x192x32/'; scaling='strong'; -% DIR = '/misc/gyacomo23_outputs/scaling/strong/3x2x768x192x32/'; scaling='strong'; +DIR = '/misc/gyacomo23_outputs/scaling/strong/3x2x768x192x32/'; scaling='strong'; % DIR = '/misc/gyacomo23_outputs/scaling/weak/Np_5x2x128x64x32/'; scaling='weak'; % DIR = '/misc/gyacomo23_outputs/scaling/weak/Ny_3x2x768x8x32/'; scaling='weak'; % DIR = '/misc/gyacomo23_outputs/scaling/weak/Nz_3x2x768x32x8/'; scaling='weak'; % DIR = '/misc/gyacomo23_outputs/scaling/weak/Nz_5x2x128x64x8/'; scaling='weak'; -% DIR = '/misc/gyacomo23_outputs/scaling/weak/Npyz_4x2x32x16x16/'; scaling='weak'; +% DIR = '/misc/gyacomo23_outputs/scaling/weak/Npyz_4x2x32x16x8/'; scaling='weak'; +% DIR = '/misc/gyacomo23_outputs/scaling/weak/Npyz_8x2x32x32x4/'; scaling='weak'; % Get a list of all items in the current directory contents = dir(DIR); @@ -76,31 +77,52 @@ Rt_ref = Rt_avg(Np_tot==1); Np_ref = 1; MaxN = max(Np_tot); Np1Nz1 = logical((Np==1).*(Nz==1)); -Np2Nz1 = logical((Np==2).*(Nz==1)); Np1Nz2 = logical((Np==1).*(Nz==2)); Np1Nz4 = logical((Np==1).*(Nz==4)); +Np2Nz1 = logical((Np==2).*(Nz==1)); +Np2Nz2 = logical((Np==2).*(Nz==2)); +Np2Nz4 = logical((Np==2).*(Nz==4)); +Np4Nz1 = logical((Np==4).*(Nz==1)); +Np4Nz2 = logical((Np==4).*(Nz==2)); +Np4Nz4 = logical((Np==4).*(Nz==4)); %% -figure; hold on; +figure; hold on; msz = 8; xlabel('Number of cores'); +clr_ = lines(3); switch scaling case 'strong' - plot(Np_tot(Np1Nz1),Rt_ref./Rt_avg(Np1Nz1),... - 'o','DisplayName','$N_{pp}=1, N_{pz}=1$'); - plot(Np_tot(Np2Nz1),Rt_ref./Rt_avg(Np2Nz1), ... - 'o','DisplayName','$N_{pp}=2, N_{pz}=1$'); - plot(Np_tot(Np1Nz2),Rt_ref./Rt_avg(Np1Nz2), ... - 'o','DisplayName','$N_{pp}=1, N_{pz}=2$'); - plot(Np_tot(Np1Nz4),Rt_ref./Rt_avg(Np1Nz4), ... - 'o','DisplayName','$N_{pp}=1, N_{pz}=4$'); - plot(1:MaxN,(1:MaxN)/Np_ref,'--k','DisplayName','Ideal'); + plot(Np_tot(Np1Nz1),Rt_ref./Rt_avg(Np1Nz1),'o', ... + 'Color',clr_(1,:),'DisplayName','$N_{pp}=1, N_{pz}=1$','MarkerSize',msz); + plot(Np_tot(Np1Nz2),Rt_ref./Rt_avg(Np1Nz2),'o',... + 'Color',clr_(2,:),'DisplayName','$N_{pp}=1, N_{pz}=2$','MarkerSize',msz); + plot(Np_tot(Np1Nz4),Rt_ref./Rt_avg(Np1Nz4),'o', ... + 'Color',clr_(3,:),'DisplayName','$N_{pp}=1, N_{pz}=4$','MarkerSize',msz); + plot(Np_tot(Np2Nz1),Rt_ref./Rt_avg(Np2Nz1),'x', ... + 'Color',clr_(1,:),'DisplayName','$N_{pp}=2, N_{pz}=1$','MarkerSize',msz); + plot(Np_tot(Np2Nz2),Rt_ref./Rt_avg(Np2Nz2),'x', ... + 'Color',clr_(2,:),'DisplayName','$N_{pp}=2, N_{pz}=2$','MarkerSize',msz); + plot(Np_tot(Np2Nz4),Rt_ref./Rt_avg(Np2Nz4),'x', ... + 'Color',clr_(3,:),'DisplayName','$N_{pp}=2, N_{pz}=4$','MarkerSize',msz); + plot(Np_tot(Np4Nz1),Rt_ref./Rt_avg(Np4Nz1),'+', ... + 'Color',clr_(1,:),'DisplayName','$N_{pp}=4, N_{pz}=1$','MarkerSize',msz); + plot(Np_tot(Np4Nz2),Rt_ref./Rt_avg(Np4Nz2),'+', ... + 'Color',clr_(2,:),'DisplayName','$N_{pp}=4, N_{pz}=2$','MarkerSize',msz); + plot(Np_tot(Np4Nz4),Rt_ref./Rt_avg(Np4Nz4),'+', ... + 'Color',clr_(3,:),'DisplayName','$N_{pp}=4, N_{pz}=4$','MarkerSize',msz); + plot(Np_tot,Rt_ref./Rt_avg,... + '.k','DisplayName','Strong scaling'); + plot([1 1000],[1 1000],'--k','DisplayName','Ideal'); + axis equal set(gca,'XScale','log'); set(gca,'YScale','log'); ylabel('Speedup'); % title('Strong scaling speedup') case 'weak' hold on; - filt = logical((Np==2).*(Ny>=1).*(Nz>=1)); - plot(Np_tot(filt),Rt_ref./Rt_avg(filt),... - 'o','DisplayName','Effective'); + % for n = [1 2 4 8 12] + filt = logical((Np>=1).*(Ny>=1).*(Nz>=1)); + plot(Np_tot(filt),Rt_ref./Rt_avg(filt),... + 'o','DisplayName','Effective'); + % end plot(1:MaxN,ones(1,MaxN),'--k','DisplayName','Ideal'); ylim([0 1]); set(gca,'XScale','log') diff --git a/wk/parameters/lin_DIIID_LM_rho95.m b/wk/parameters/lin_DIIID_LM_rho95.m index 2bd5ba9b2f383bfb8a885ecb30666f6b79194743..17c6cbf8552518eb05c6cbd2b62d9671e59d0730 100644 --- a/wk/parameters/lin_DIIID_LM_rho95.m +++ b/wk/parameters/lin_DIIID_LM_rho95.m @@ -96,4 +96,5 @@ BCKGD0 = 0.0e-5; % Initial background k_gB = 1.0; % Magnetic gradient strength k_cB = 1.0; % Magnetic curvature strength COLL_KCUT = 1; % Cutoff for collision operator -ADIAB_I = 0; % adiabatic ion model \ No newline at end of file +ADIAB_I = 0; % adiabatic ion model +EXBRATE = 0; \ No newline at end of file diff --git a/wk/parameters/lin_DIIID_LM_rho95_HEL.m b/wk/parameters/lin_DIIID_LM_rho95_HEL.m new file mode 100644 index 0000000000000000000000000000000000000000..cbd940901c56c65a2f11ccf80320324342097e03 --- /dev/null +++ b/wk/parameters/lin_DIIID_LM_rho95_HEL.m @@ -0,0 +1,100 @@ +%% Reference values +% See Neiser et al. 2019 Gyrokinetic GENE simulations of DIII-D near-edge L-mode plasmas +%% Set simulation parameters +SIMID = 'lin_DIIID_HM_rho90'; % Name of the simulation +%% Set up physical parameters +CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss +NU = 50; %(0.00235 in GENE) +TAU = 0.001; % i/e temperature ratio +K_Ne = 0; % ele Density ''' +K_Te = 0; % ele Temperature ''' +K_Ni = 0; % ion Density gradient drive +K_Ti = 5.2256/2/TAU; % ion Temperature ''' +SIGMA_E = 0.0233380/sqrt(2); % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) +NA = 1; % number of kinetic species +ADIAB_E = (NA==1); % adiabatic electron model +BETA = 0; % electron plasma beta in prct +MHD_PD = 1; +%% Set up grid parameters +P = 2; +J = P/2;%P/2; +PMAX = P; % Hermite basis size +JMAX = J; % Laguerre basis size +NX = 8; % real space x-gridpoints +NY = 2; % real space y-gridpoints +LX = 2*pi/0.1; % Size of the squared frequency domain in x direction +LY = 2*pi/0.3; % Size of the squared frequency domain in y direction +NZ = 32; % number of perpendicular planes (parallel grid) +SG = 0; % Staggered z grids option +NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) + +%% GEOMETRY +% GEOMETRY= 's-alpha'; +GEOMETRY= 'miller'; +Q0 = 4.8; % safety factor +SHEAR = 2.55; % magnetic shear +EPS = 0.3; % inverse aspect ratio +KAPPA = 1.57; % elongation +S_KAPPA = 0.48; +DELTA = 0.2; % triangularity +S_DELTA = 0.1; +ZETA = 0.0; % squareness +S_ZETA = 0.0; +PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected') +SHIFT_Y = 0.0; % Shift in the periodic BC in z +NPOL = 1; % Number of poloidal turns +PB_PHASE = 0; +%% TIME PARAMETERS +TMAX = 15; % Maximal time unit +DT = 1e-2; % Time step +DTSAVE0D = 0.5; % Sampling time for 0D arrays +DTSAVE2D = -1; % Sampling time for 2D arrays +DTSAVE3D = 0.5; % Sampling time for 3D arrays +DTSAVE5D = 100; % Sampling time for 5D arrays +JOB2LOAD = -1; % Start a new simulation serie + +%% OPTIONS +LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) +CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) +GKCO = 1; % Gyrokinetic operator +ABCO = 1; % INTERSPECIES collisions +INIT_ZF = 0; % Initialize zero-field quantities +HRCY_CLOS = 'truncation'; % Closure model for higher order moments +DMAX = -1; +NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments +NMAX = 0; +KERN = 0; % Kernel model (0 : GK) +INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom +NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5) + +%% OUTPUTS +W_DOUBLE = 1; % Output flag for double moments +W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy) +W_HF = 1; % Output flag for high-frequency potential energy +W_PHI = 1; % Output flag for potential +W_NA00 = 1; % Output flag for nalpha00 (density of species alpha) +W_DENS = 1; % Output flag for total density +W_TEMP = 1; % Output flag for temperature +W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha) +W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% UNUSED PARAMETERS +% These parameters are usually not to play with in linear runs +MU = 0.0; % Hyperdiffusivity coefficient +MU_X = MU; % Hyperdiffusivity coefficient in x direction +MU_Y = MU; % Hyperdiffusivity coefficient in y direction +N_HD = 4; % Degree of spatial-hyperdiffusivity +MU_Z = 5.0; % Hyperdiffusivity coefficient in z direction +HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model +MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite +MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre +LAMBDAD = 0.0; % Lambda Debye +NOISE0 = 1.0e-5; % Initial noise amplitude +BCKGD0 = 0.0e-5; % Initial background +k_gB = 1.0; % Magnetic gradient strength +k_cB = 1.0; % Magnetic curvature strength +COLL_KCUT = 1; % Cutoff for collision operator +ADIAB_I = 0; % adiabatic ion model +EXBRATE = 0; \ No newline at end of file diff --git a/wk/parameters/lin_ITG.m b/wk/parameters/lin_ITG.m index 07d85ebc3aa01595ebd7e924280e8eaf25d875c9..d1c06b664d3094daa98e7e1c8e4c327e342cf16f 100644 --- a/wk/parameters/lin_ITG.m +++ b/wk/parameters/lin_ITG.m @@ -15,10 +15,8 @@ ADIAB_E = (NA==1); % adiabatic electron model BETA = 0.0; % electron plasma beta EXBRATE = 0.0; % Background ExB shear flow %% Set up grid parameters -P = 4; -J = 2;%P/2; -PMAX = P; % Hermite basis size -JMAX = J; % Laguerre basis size +PMAX = 4; % Hermite basis size +JMAX = 2; % Laguerre basis size NX = 8; % real space x-gridpoints NY = 12; % real space y-gridpoints LX = 2*pi/0.05; % Size of the squared frequency domain in x direction diff --git a/wk/parameters/lin_Ivanov.m b/wk/parameters/lin_Ivanov.m index 80dccb8e179b3b7fae17dc3faa49e3eef211c83d..a578f0bc8ec6a04a2ea32e0cd771d7dd45ab8f1b 100644 --- a/wk/parameters/lin_Ivanov.m +++ b/wk/parameters/lin_Ivanov.m @@ -3,12 +3,12 @@ SIMID = 'lin_Ivanov'; % Name of the simulation %% Set up physical parameters CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss -TAU = 0.001; % e/i temperature ratio -NU = 1.0*3/2/TAU/4; % Collision frequency +TAU = 1e-3; % e/i temperature ratio +NU = 0.1*3/8/TAU; % Collision frequency K_Ne = 0*2.22; % ele Density K_Te = 0*6.96; % ele Temperature K_Ni = 0*2.22; % ion Density gradient drive -K_Ti = 1.0*2/TAU+0*5*TAU; % ion Temperature +K_Ti = 1*2/TAU; % ion Temperature SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) NA = 1; % number of kinetic species ADIAB_E = (NA==1); % adiabatic electron model @@ -96,7 +96,10 @@ BCKGD0 = 0.0; % Initial background k_gB = 1.0; % Magnetic gradient strength k_cB = 1.0; % Magnetic curvature strength COLL_KCUT = 1.0; % Cutoff for collision operator -PB_PHASE = 0.0; +S_KAPPA = 0.0; +S_DELTA = 0.0; +S_ZETA = 0.0; +PB_PHASE= 0; ADIAB_I = 0; MHD_PD = 0; EXBRATE = 0; \ No newline at end of file diff --git a/wk/parameters/lin_RHT.m b/wk/parameters/lin_RHT.m index 72137f75d6cb2f219d7ce3bf3bd81ecc77aea40e..0f39498d342d28c11b63c43d64e70132d1f07e0a 100644 --- a/wk/parameters/lin_RHT.m +++ b/wk/parameters/lin_RHT.m @@ -15,15 +15,15 @@ ADIAB_E = (NA==1); % adiabatic electron model BETA = 0.0; % electron plasma beta EXBRATE = 0.0; % Background ExB shear flow %% Set up grid parameters -P = 4; -J = 2;%P/2; +P = 256; +J = 16;%P/2; PMAX = P; % Hermite basis size JMAX = J; % Laguerre basis size NX = 2; % real space x-gridpoints NY = 2; % real space y-gridpoints LX = 2*pi/0.05; % Size of the squared frequency domain in x direction LY = 2*pi/0.01; % Size of the squared frequency domain in y direction -NZ = 8; % number of perpendicular planes (parallel grid) +NZ = 24; % number of perpendicular planes (parallel grid) SG = 0; % Staggered z grids option NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) @@ -42,7 +42,7 @@ NPOL = 1; % Number of poloidal turns %% TIME PARAMETERS TMAX = 50; % Maximal time unit -DT = 1e-3; % Time step +DT = 5e-3; % Time step DTSAVE0D = 1; % Sampling time for 0D arrays DTSAVE2D = -1; % Sampling time for 2D arrays DTSAVE3D = 0.1; % Sampling time for 3D arrays diff --git a/wk/particle_trajectory.m b/wk/particle_trajectory.m new file mode 100644 index 0000000000000000000000000000000000000000..62da611a8b43f50d5b00184918e8d22416b0b2ed --- /dev/null +++ b/wk/particle_trajectory.m @@ -0,0 +1,36 @@ +R = 1; +a = 0.05; +d = 0.02; +Wc= 80; + + +figure + +% Magnetic field line +Bx = @(R,th) R*cos(th); +By = @(R,th) R*sin(th); +Bz = @(R,th) 0*(th); + +thB = linspace(0,2*pi,128); +plot3(Bx(R,thB),By(R,thB),Bz(R,thB)); hold on; + +% Particle trochoid +Tx = @(R,th) a*cos(th).*cos(Wc*th); +Ty = @(R,th) a*sin(th).*cos(Wc*th); +Tz = @(R,th) a*sin(Wc*th)+d*th; + +thP = linspace(0,2.5*pi,2048); +Xp = Bx(R,thP)+Tx(R,thP); +Yp = By(R,thP)+Ty(R,thP); +Zp = Bz(R,thP)+Tz(R,thP); +plot3(Xp,Yp,Zp) +plot3(Xp(1),Yp(1),Zp(1),'xk','MarkerSize',10) +plot3(Xp(end),Yp(end),Zp(end),'ok','MarkerSize',8,'MarkerFaceColor','k') + + +% finish the plot +axis equal +xlim(1.2*R*[-1 1]) +ylim(1.2*R*[-1 1]) +zlim(0.25*[-1 1]) +axis off \ No newline at end of file diff --git a/wk/plot_params_vs_rho.m b/wk/plot_params_vs_rho.m index 0b3d64860fbba428210049054fa658e9b577785f..0f12e391a80a4bc72a373a40bc098cdb4bec93ae 100644 --- a/wk/plot_params_vs_rho.m +++ b/wk/plot_params_vs_rho.m @@ -3,6 +3,7 @@ function [ ] = plot_params_vs_rho(geom,prof_folder,rho,rho_min,rho_max,Lref,mref if FROMPROFILE % profiles = read_DIIID_profile([prof_folder,'/profiles.txt']); [params,profiles] = get_param_from_profiles(prof_folder,rho,Lref,mref,Bref,FROMPROFILE); + rho_a = profiles.ne.x; m_e = 9.11e-31; sigma = sqrt(m_e/mref); TAU_a = profiles.ti.y./profiles.te.y; @@ -42,12 +43,15 @@ if FROMPROFILE else rho_a = linspace(rho_min,rho_max,100); - NU_a = zeros(size(rho_a)); - BETA_a = zeros(size(rho_a)); - K_Ne_a = zeros(size(rho_a)); - K_Ti_a = zeros(size(rho_a)); - K_Te_a = zeros(size(rho_a)); - + NU_a = zeros(size(rho_a)); + BETA_a = zeros(size(rho_a)); + K_Ne_a = zeros(size(rho_a)); + K_Ti_a = zeros(size(rho_a)); + K_Te_a = zeros(size(rho_a)); + Kappa_a = zeros(size(rho_a)); + Delta_a = zeros(size(rho_a)); + Zeta_a = zeros(size(rho_a)); + geom_a = cell(size(rho_a)); for i = 1:numel(rho_a) rho_ = rho_a(i); [params,profiles] = get_param_from_profiles(prof_folder,rho_,Lref,mref,Bref,FROMPROFILE); @@ -56,6 +60,12 @@ else K_Ne_a(i) = params.K_Ne; K_Ti_a(i) = params.K_Ti; K_Te_a(i) = params.K_Te; + + geom_ = get_miller_GENE_py(prof_folder,rho_); + Kappa_a(i) = geom_.kappa; + Delta_a(i) = geom_.delta; + Zeta_a(i) = geom_.zeta; + geom_a{i} = geom_; end figure subplot(231) @@ -88,8 +98,15 @@ else grid on end subplot(133) - [R,Z] = plot_miller(geom,rho,32,0); - plot(R,Z,'-b'); + [R,Z] = plot_miller(geom,rho,128,0); + plot(R,Z,'-b'); hold on + % rhos = linspace(0.09,0.99,6); + % for i = 1:numel(rhos) + % rho_ = rhos(i); + % geom_ = get_miller_GENE_py(prof_folder,rho_); + % [R,Z] = plot_miller(geom_,rho_,128,0); + % plot(R,Z,'-k'); + % end xlabel('R [m]'); ylabel('Z [m]'); axis tight diff --git a/wk/plot_spectrum.m b/wk/plot_spectrum.m new file mode 100644 index 0000000000000000000000000000000000000000..264f3a72b3acb75aaf3e96e22353b17192cf3b6f --- /dev/null +++ b/wk/plot_spectrum.m @@ -0,0 +1,45 @@ +function [FIGURE] = plot_spectrum(data,options) + + +options.INTERP = 0; +options.POLARPLOT = 0; +options.AXISEQUAL = 1; +options.LOGSCALE = 0; +options.CLIMAUTO = 1; +options.PLAN = 'kxky'; options.COMP = data.grids.Nz/2+1; +options.RESOLUTION = 256; +options.BWR = 0; % bluewhitered plot or gray +options.COMP = 'avg'; +[toplot] = process_field(data,options); + +% Time averaging +fkykx = squeeze(mean(abs(toplot.FIELD(:,:,:)),3)); + +FNAME = toplot.FILENAME; + + +kx = toplot.X(1,:); +Nkx = numel(kx); +kx = toplot.X(1,Nkx/2:end); +fkx = squeeze(fkykx(1,Nkx/2:end)); +ky = toplot.Y(:,data.grids.Nkx/2+1); +fky = squeeze(fkykx(:,data.grids.Nkx/2+1)); + +FIGURE.fig = figure; +subplot(121) + scale = max(fkx)*options.NORMALIZE + 1*(1-options.NORMALIZE); + plot(kx,fkx./scale) + xlabel(toplot.XNAME) + ylabel(['$|',toplot.FIELDNAME,'|$']) + + +subplot(122) + scale = max(fky)*options.NORMALIZE + 1*(1-options.NORMALIZE); + plot(ky,fky./scale) + xlabel(toplot.YNAME) + ylabel(['$|',toplot.FIELDNAME,'|$']) + + +FIGURE.FIGNAME = [FNAME,'_spectrum']; + +end \ No newline at end of file diff --git a/wk/plot_toroidal_fluxtube_geom.m b/wk/plot_toroidal_fluxtube_geom.m index b6aa8e4e1a8a7c0a381a08da32cd94842ca0fd65..f38037a99451e2dc7b3820340302ef82fd143510 100644 --- a/wk/plot_toroidal_fluxtube_geom.m +++ b/wk/plot_toroidal_fluxtube_geom.m @@ -1,3 +1,8 @@ +%% This script is made for illustrating different magnetic geometry +% it is not bugfree but work pretty well for illustration purpose of the +% magnetic flux surface of a tokamak and for the particle trajectory in it. +% ! The particle trajectory is not correct (does not work for Z-pinch for +% ! example) but is good enough for illustation. gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory addpath(genpath([gyacomodir,'matlab'])) % ... add addpath(genpath([gyacomodir,'matlab/plot'])) % ... add @@ -8,27 +13,31 @@ default_plots_options Nx = 128; Ny = 128; -Nturns = 1; -Nz = 128; -x = linspace(-60,60,Nx)*0.001; +Nz = 64; +Nturns = 1.0; +lr = 0.05; %larmor radius +Wc = 750; %~cyclotron freq +vd = 0.25; %~drift velocty +Gx = linspace(-60,60,Nx)*0.001; y = linspace(-60,60,Ny)*0.001; FIELD = ones(Nx,Ny,Nz); z = linspace(-Nturns*pi,Nturns*pi,Nz); -N_field_lines = 2; +N_field_lines = 10; N_magn_flux_surf = 1; openangle = pi/3; -phi0 = openangle/2; phi1= 2*pi-openangle/2; +% phi0 = openangle/2; phi1= 2*pi-openangle/2; +phi0 = 0; phi1= 2*pi-1.3; PLT_FTUBE = 0; PLT_BASIS = 0; PLT_DATA = 0; R0 = 1; % Torus major radius Z0 = 0; xpoint = 0; -select = 3; +select = 1; switch select case 1 - % CBC - rho = 0.3; drho = 0.1;% Torus minor radius + % Z-pinch + rho = 0.50; drho = 0.1;% Torus minor radius eps = 0.3; q0 = 1.4; % Flux tube safety factor shear = 0.8; @@ -78,36 +87,20 @@ switch select zeta = 0.1; s_zeta = 0.0; case 5 - % Fake x-point DIII-D - N_magn_flux_surf = 2; - phi0 = 0; phi1= pi/48; - N_field_lines = 0; - rho = 0.40; drho = 0.1;% Torus minor radius - eps = 0.3; - q0 = 4.8; % Flux tube safety factor - shear = 2.5; - kappa = 1.0; - s_kappa= 0.0; - delta = 0.0; - s_delta= 0.0; - zeta = 0.0; - s_zeta = 0.0; - xpoint = 0.5; - case 6 % play with DIII-D edge rho = 0.90; drho = 0.1;% Torus minor radius eps = 0.3; - q0 = 1.0; % Flux tube safety factor - Nturns = max(1,q0); + q0 = 1.8; % Flux tube safety factor + % Nturns = max(1,q0); shear = 2.5; - kappa = 0.5; + kappa = 1.5; s_kappa= 1.0; - delta = 0.0; + delta = 0.3; s_delta= 0.0; - zeta =-1.0; + zeta = 0.0; s_zeta = 0.0; end - +% Nturns = 1/q0; Nptor = 64; % Toroidal angle for the flux surf % phi = linspace(0+openangle/2, 2*pi-openangle/2, Nptor); @@ -115,20 +108,23 @@ phi = linspace(phi0,phi1, Nptor); theta = linspace(-pi, pi, Nptor) ; % Poloidal angle [p, t] = meshgrid(phi, theta); Tgeom = 1; - DIMENSIONS = [600 600 1200 600]; rho = rho*eps; drho = drho*eps; +iota = 1/q0; % field line coordinates Xfl = @(z) (R0+rho*cos(z+asin(delta*sin(z)))).*cos(q0*z); Yfl = @(z) (R0+rho*cos(z+asin(delta*sin(z)))).*sin(q0*z); Zfl = @(z) Z0+kappa*rho*sin(z+zeta*sin(2*z)+xpoint*cos(2*z)); Rvec= @(z) [Xfl(z); Yfl(z); Zfl(z)]; % xvec shearless -xX = @(z) (Xfl(z)-R0*cos(q0*z))./sqrt((Xfl(z)-R0*cos(q0*z)).^2+(Yfl(z)-R0*sin(q0*z)).^2+Zfl(z).^2); -xY = @(z) (Yfl(z)-R0*sin(q0*z))./sqrt((Xfl(z)-R0*cos(q0*z)).^2+(Yfl(z)-R0*sin(q0*z)).^2+Zfl(z).^2); -xZ = @(z) Zfl(z)./sqrt((Xfl(z)-R0*cos(q0*z)).^2+(Yfl(z)-R0*sin(q0*z)).^2+Zfl(z).^2); +% xX = @(z) (Xfl(z)-R0*cos(q0*z))./sqrt((Xfl(z)-R0*cos(q0*z)).^2+(Yfl(z)-R0*sin(q0*z)).^2+Zfl(z).^2); +% xY = @(z) (Yfl(z)-R0*sin(q0*z))./sqrt((Xfl(z)-R0*cos(q0*z)).^2+(Yfl(z)-R0*sin(q0*z)).^2+Zfl(z).^2); +% xZ = @(z) Zfl(z)./sqrt((Xfl(z)-R0*cos(q0*z)).^2+(Yfl(z)-R0*sin(q0*z)).^2+Zfl(z).^2); +xX = @(z) cos(z+asin(delta*sin(z))).*cos(q0*z)./sqrt(cos(z+asin(delta*sin(z))).^2+kappa^2*sin(z+zeta*sin(2*z)).^2); +xY = @(z) cos(z+asin(delta*sin(z))).*sin(q0*z)./sqrt(cos(z+asin(delta*sin(z))).^2+kappa^2*sin(z+zeta*sin(2*z)).^2); +xZ = @(z) kappa*sin(z+zeta*sin(2*z))./sqrt(cos(z+asin(delta*sin(z))).^2+kappa^2*sin(z+zeta*sin(2*z)).^2); xvec= @(z) [xX(z); xY(z); xZ(z)]; -% bvec +% bvec shearless bX = @(z) Tgeom*(rho*cos(z).*cos(q0*z) - q0*Yfl(z))./sqrt(Tgeom*(rho*cos(z).*cos(q0*z) - q0*Yfl(z)).^2+(rho*cos(z).*sin(q0*z) + q0*Xfl(z)).^2+(rho*cos(z)).^2); bY = @(z) (rho*cos(z).*sin(q0*z) + q0*Xfl(z))./sqrt(Tgeom*(rho*cos(z).*cos(q0*z) - q0*Yfl(z)).^2+(rho*cos(z).*sin(q0*z) + q0*Xfl(z)).^2+(rho*cos(z)).^2); bZ = @(z) rho*cos(z)./sqrt(Tgeom*(rho*cos(z).*cos(q0*z) - q0*Yfl(z)).^2+(rho*cos(z).*sin(q0*z) + q0*Xfl(z)).^2+(rho*cos(z)).^2); @@ -141,13 +137,6 @@ yvec= @(z) [yX(z); yY(z); yZ(z)]; scale = 0.10; -% Plot plane result -OPTIONS.POLARPLOT = 0; -OPTIONS.PLAN = 'xy'; -rhostar = 0.1; -[X,Y] = meshgrid(x,y); -max_ = 0; - figure; set(gcf, 'Position', DIMENSIONS) %plot magnetic geometry @@ -171,7 +160,8 @@ figure; set(gcf, 'Position', DIMENSIONS) end %plot field lines if N_field_lines > 0 - theta = linspace(-Nturns*pi, Nturns*pi, Nturns*256) ; % Poloidal angle + % theta = linspace(-Nturns*pi, Nturns*pi, 1024) ; % Poloidal angle + theta = linspace(0, Nturns*2*pi, 1024) ; % Poloidal angle vecfl = Rvec(theta); plot3(vecfl(1,:),vecfl(2,:),vecfl(3,:),'-b'); hold on; % Multiple field lines @@ -179,17 +169,22 @@ figure; set(gcf, 'Position', DIMENSIONS) for ifl = 1:N_field_lines-1 % rotation for multiple field lines t = q0*pi*ifl; - x_ = [x_ cos(t)*vecfl(1,:) - sin(t)*vecfl(2,:)]; - y_ = [y_ sin(t)*vecfl(1,:) + cos(t)*vecfl(2,:)]; - z_ = [z_ vecfl(3,:)]; + % x_ = [x_ cos(t)*vecfl(1,:) - sin(t)*vecfl(2,:)]; + % y_ = [y_ sin(t)*vecfl(1,:) + cos(t)*vecfl(2,:)]; + % z_ = [z_ vecfl(3,:)]; + x_ = cos(t)*vecfl(1,:) - sin(t)*vecfl(2,:); + y_ = sin(t)*vecfl(1,:) + cos(t)*vecfl(2,:); + z_ = vecfl(3,:); + % plot3(x_,y_,z_,'-','color',[1.0 0.6 0.6]*0.8); hold on; + plot3(x_,y_,z_,'-','color',[0.0 0.5 1.0]*0.8,'LineWidth',1); hold on; end - plot3(x_,y_,z_,'-','color',[1.0 0.6 0.6]*0.8); hold on; + %plot3(x_,y_,z_,'-','color',[1.0 0.6 0.6]*0.8); hold on; end %plot fluxe tube if PLT_FTUBE theta = linspace(-Nturns*pi, Nturns*pi, 64) ; % Poloidal angle %store the shifts in an order (top left to bottom right) - s_x = rhostar*[x(1) x(end) x(1) x(end)]; + s_x = rhostar*[Gx(1) Gx(end) Gx(1) Gx(end)]; s_y = rhostar*[y(1) y(1) y(end) y(end)]; for i_ = 1:4 vx_ = Xfl(theta) + s_x(i_)*xX(theta) + s_y(i_)*yX(theta); @@ -206,19 +201,49 @@ figure; set(gcf, 'Position', DIMENSIONS) quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*yX(theta),scale*yY(theta),scale*yZ(theta),0,'g'); quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*bX(theta),scale*bY(theta),scale*bZ(theta),0,'b'); end - xlabel('X');ylabel('Y');zlabel('Z'); - %Plot time dependent data - if PLT_DATA - for iz = 1:Nz - z_ = z(iz); - Xp = Xfl(z_) + X*xX(z_) + Y*yX(z_); - Yp = Yfl(z_) + X*xY(z_) + Y*yY(z_); - Zp = Zfl(z_) + X*xZ(z_) + Y*yZ(z_); - s=surface(Xp,Yp,Zp,FIELD(:,:,iz)/max(max(max(abs(FIELD))))); - s.EdgeColor = 'none'; - colormap(bluewhitered); + % Particle trajectory + if 1 + % Particle gyration + Gx = @(z) lr*cos(z).*cos(Wc*z); + Gy = @(z) lr*sin(z).*cos(Wc*z); + Gz = @(z) lr*sin(Wc*z); + % This is the "real" helicoidal trajectory + theta = linspace(0, Nturns*2*pi, 10000) ; % Poloidal angle + dtheta = theta(2)-theta(1); + dX = 0*theta; dY = dX; dZ = dX; % init drifts + for i = 2:numel(theta) + th_old = theta(i-1); + th_ = theta(i); + % find magnetic displacement vector + bx_= Xfl(th_)-Xfl(th_old); + by_= Yfl(th_)-Yfl(th_old); + bz_= Zfl(th_)-Zfl(th_old); + % find binormal + yx_= xY(th_).*bz_ - xZ(th_).*by_; + yy_= xZ(th_).*bx_ - xX(th_).*bz_; + yz_= xX(th_).*by_ - xY(th_).*bx_; + % normalize + yx_=yx_/sqrt(yx_^2+yy_^2+yz_^2); + yy_=yy_/sqrt(yx_^2+yy_^2+yz_^2); + yz_=yz_/sqrt(yx_^2+yy_^2+yz_^2); + % add drift + dX(i) = dX(i-1) - vd*yx_*dtheta; + dY(i) = dY(i-1) - vd*yy_*dtheta; + dZ(i) = dZ(i-1) - vd*yz_*dtheta; end + % add gyration + dX = dX + Gx(theta); + dY = dY + Gy(theta); + dZ = dZ + Gz(theta); + % add fieldline + Tx = Xfl(theta)+dX; + Ty = Yfl(theta)+dY; + Tz = Zfl(theta)+dZ; + plot3(Tx,Ty,Tz,'-r','LineWidth',1.2); hold on; + plot3(Tx(1),Ty(1),Tz(1),'xk','MarkerSize',10); hold on; + plot3(Tx(end),Ty(end),Tz(end),'ok','MarkerSize',8,'MarkerFaceColor','k'); hold on; end + xlabel('X');ylabel('Y');zlabel('Z'); % axis equal view([1,-2,1])