diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m index 6f5e640..c468731 100644 --- a/matlab/plot/plot_radial_transport_and_spacetime.m +++ b/matlab/plot/plot_radial_transport_and_spacetime.m @@ -1,143 +1,145 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS) %Compute steady radial transport tend = OPTIONS.TAVG_1; tstart = OPTIONS.TAVG_0; [~,its0D] = min(abs(DATA.Ts0D-tstart)); [~,ite0D] = min(abs(DATA.Ts0D-tend)); [~,its3D] = min(abs(DATA.Ts3D-tstart)); [~,ite3D] = min(abs(DATA.Ts3D-tend)); SCALE = DATA.scale;%(1/DATA.Nx/DATA.Ny)^2; Gx_infty_avg = mean(DATA.PGAMMA_RI(its0D:ite0D))*SCALE; Gx_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE; Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE; Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE; disp(['G_x=',sprintf('%2.2e',Gx_infty_avg),'+-',sprintf('%2.2e',Gx_infty_std)]); % disp(['Q_x=',sprintf('%2.2e',Qx_infty_avg),'+-',sprintf('%2.2e',Qx_infty_std)]); f_avg_z = squeeze(mean(DATA.PHI(:,:,:,:),3)); [~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3))); ikzf = min([ikzf,DATA.Nky]); Ns3D = numel(DATA.Ts3D); [KX, KY] = meshgrid(DATA.kx, DATA.ky); %% error estimation DT_ = (tend-tstart)/OPTIONS.NCUT; Qx_ee = zeros(1,OPTIONS.NCUT); for i = 1:OPTIONS.NCUT [~,its_] = min(abs(DATA.Ts0D - (tstart+(i-1)*DT_))); [~,ite_] = min(abs(DATA.Ts0D - (tstart+ i *DT_))); Qx_ee(i) = mean(DATA.HFLUX_X(its_:ite_))*SCALE; end Qx_avg = mean(Qx_ee); Qx_err = std(Qx_ee); % disp(['Q_avg=',sprintf('%2.2e',Qx_avg),'+-',sprintf('%2.2e',Qx_err)]); %% computations % Compute Gamma from ifft matlab Gx = zeros(DATA.Ny,DATA.Nx,numel(DATA.Ts3D)); % for it = 1:numel(DATA.Ts3D) % for iz = 1:DATA.Nz % Gx(:,:,it) = Gx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))... % .*ifourier_GENE(DATA.DENS_I(:,:,iz,it)); % end % Gx(:,:,it) = Gx(:,:,it)/DATA.Nz; % end Gx_t_mtlb = squeeze(mean(mean(Gx,1),2)); % Compute Heat flux from ifft matlab Qx = zeros(DATA.Ny,DATA.Nx,numel(DATA.Ts3D)); % for it = 1:numel(DATA.Ts3D) % for iz = 1:DATA.Nz % Qx(:,:,it) = Qx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))... % .*ifourier_GENE(DATA.TEMP_I(:,:,iz,it)); % end % Qx(:,:,it) = Qx(:,:,it)/DATA.Nz; % end Qx_t_mtlb = squeeze(mean(mean(Qx,1),2)); % zonal vs nonzonal energies for phi(t) E_Zmode_SK = zeros(1,Ns3D); E_NZmode_SK = zeros(1,Ns3D); for it = 1:numel(DATA.Ts3D) E_Zmode_SK(it) = squeeze(DATA.ky(ikzf).^2.*abs(squeeze(f_avg_z(ikzf,1,it))).^2); E_NZmode_SK(it) = squeeze(sum(sum(((1+KX.^2+KY.^2).*abs(squeeze(f_avg_z(:,:,it))).^2.*(KY~=0))))); end %% Figure mvm = @(x) movmean(x,OPTIONS.NMVA); FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1000, 600]) subplot(311) % yyaxis left plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on; % plot(mvm(DATA.Ts3D),mvm(Gx_t_mtlb),'DisplayName','matlab comp.'); hold on; % plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Gx_infty_avg, '-k',... % 'DisplayName',['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$\pm$',num2str(Gx_infty_std)]); % ylabel('$\Gamma_x$') % ylim([0,5*abs(Gx_infty_avg)]); % xlim([DATA.Ts0D(1),DATA.Ts0D(end)]); % yyaxis right plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on; % plot(mvm(DATA.Ts3D),mvm(Qx_t_mtlb),'DisplayName','matlab comp.'); hold on; ylabel('Transport') if(~isnan(Qx_infty_avg)) plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',... 'DisplayName',['$Q_{avg}=',sprintf('%2.2f',Qx_avg),'\pm',sprintf('%2.2f',Qx_err),'$']); legend('show'); ylim([0,5*abs(Qx_infty_avg)]); else plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Gx_infty_avg, '-k',... 'DisplayName',['$\Gamma_{avg}=',sprintf('%2.2f',Gx_infty_avg),'\pm',sprintf('%2.2f',Gx_infty_std),'$']); legend('show'); - ylim([0,5*abs(Gx_infty_avg)]); + try ylim([0,5*abs(Gx_infty_avg)]); + catch + end end xlim([DATA.Ts0D(1),DATA.Ts0D(end)]); 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 % computation Ns3D = numel(DATA.Ts3D); [KX, KY] = meshgrid(DATA.kx, DATA.ky); plt = @(x) mean(x(:,:,:),1); kycut = max(DATA.ky); kxcut = max(DATA.kx); LP = (abs(KY) gather phi..' CALL gather_xyz(phi(ikys:ikye,ikxs:ikxe,izs:ize), field_to_check) IF(my_id.EQ. 0) THEN WRITE(check_filename,'(a16)') 'check_phi.out' OPEN(fid_check, file=check_filename, form='formatted') WRITE(*,*) 'Check file found -> output phi ..' WRITE(fid_check,*) Nky, Nkx, Nz DO iky = 1,Nky; DO ikx = 1, Nkx; DO iz = 1,Nz WRITE(fid_check,*) real(field_to_check(iky,ikx,iz)), ',' , imag(field_to_check(iky,ikx,iz)) ENDDO; ENDDO; ENDDO CLOSE(fid_check) WRITE(*,*) 'Check file found -> done.' ! delete the check_phi flagfile OPEN(fid_check, file='check_phi') CLOSE(fid_check, status='delete') ENDIF ENDIF END SUBROUTINE spit_snapshot_check diff --git a/src/grid_mod.F90 b/src/grid_mod.F90 index 06e4710..191b319 100644 --- a/src/grid_mod.F90 +++ b/src/grid_mod.F90 @@ -1,629 +1,630 @@ MODULE grid ! Grid module for spatial discretization USE prec_const USE basic IMPLICIT NONE PRIVATE ! GRID Namelist INTEGER, PUBLIC, PROTECTED :: pmaxe = 1 ! The maximal electron Hermite-moment computed INTEGER, PUBLIC, PROTECTED :: jmaxe = 1 ! The maximal electron Laguerre-moment computed INTEGER, PUBLIC, PROTECTED :: pmaxi = 1 ! The maximal ion Hermite-moment computed INTEGER, PUBLIC, PROTECTED :: jmaxi = 1 ! The maximal ion Laguerre-moment computed INTEGER, PUBLIC, PROTECTED :: maxj = 1 ! The maximal Laguerre-moment INTEGER, PUBLIC, PROTECTED :: dmaxe = 1 ! The maximal full GF set of e-moments v^dmax INTEGER, PUBLIC, PROTECTED :: dmaxi = 1 ! The maximal full GF set of i-moments v^dmax INTEGER, PUBLIC, PROTECTED :: Nx = 16 ! Number of total internal grid points in x REAL(dp), PUBLIC, PROTECTED :: Lx = 1._dp ! horizontal length of the spatial box INTEGER, PUBLIC, PROTECTED :: Nexc = 1 ! factor to increase Lx when shear>0 (Lx = Nexc/kymin/shear) INTEGER, PUBLIC, PROTECTED :: Ny = 16 ! Number of total internal grid points in y REAL(dp), PUBLIC, PROTECTED :: Ly = 1._dp ! vertical length of the spatial box INTEGER, PUBLIC, PROTECTED :: Nz = 1 ! Number of total perpendicular planes INTEGER, PUBLIC, PROTECTED :: Npol = 1 ! number of poloidal turns INTEGER, PUBLIC, PROTECTED :: Odz = 4 ! order of z interp and derivative schemes INTEGER, PUBLIC, PROTECTED :: Nkx = 8 ! Number of total internal grid points in kx REAL(dp), PUBLIC, PROTECTED :: Lkx = 1._dp ! horizontal length of the fourier box INTEGER, PUBLIC, PROTECTED :: Nky = 16 ! Number of total internal grid points in ky REAL(dp), PUBLIC, PROTECTED :: Lky = 1._dp ! vertical length of the fourier box REAL(dp), PUBLIC, PROTECTED :: kpar = 0_dp ! parallel wave vector component ! For Orszag filter REAL(dp), PUBLIC, PROTECTED :: two_third_kxmax REAL(dp), PUBLIC, PROTECTED :: two_third_kymax REAL(dp), PUBLIC :: two_third_kpmax ! 1D Antialiasing arrays (2/3 rule) REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_x REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_y ! Grids containing position in physical space REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: xarray REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: yarray ! Local and global z grids, 2D since it has to store odd and even grids REAL(dp), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: zarray REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: zarray_full ! local z weights for computing simpson rule INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: zweights_SR REAL(dp), PUBLIC, PROTECTED :: deltax, deltay, deltaz, inv_deltaz REAL(dp), PUBLIC, PROTECTED :: diff_kx_coeff, diff_ky_coeff, diff_dz_coeff INTEGER, PUBLIC, PROTECTED :: ixs, ixe, iys, iye, izs, ize INTEGER, PUBLIC, PROTECTED :: izgs, izge ! ghosts LOGICAL, PUBLIC, PROTECTED :: SG = .true.! shifted grid flag INTEGER, PUBLIC :: ir,iz ! counters ! Data about parallel distribution for ky.kx integer(C_INTPTR_T), PUBLIC :: local_nkx, local_nky integer(C_INTPTR_T), PUBLIC :: local_nkx_offset, local_nky_offset INTEGER, PUBLIC :: local_nkp INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: counts_nkx, counts_nky INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: displs_nkx, displs_nky ! "" for p INTEGER, PUBLIC :: local_np_e, local_np_i INTEGER, PUBLIC :: total_np_e, total_np_i, Np_e, Np_i integer(C_INTPTR_T), PUBLIC :: local_np_e_offset, local_np_i_offset INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: rcv_p_e, rcv_p_i INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: dsp_p_e, dsp_p_i ! "" for z INTEGER, PUBLIC :: local_nz INTEGER, PUBLIC :: total_nz integer(C_INTPTR_T), PUBLIC :: local_nz_offset INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: counts_nz INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: displs_nz ! "" for j (not parallelized) INTEGER, PUBLIC :: local_nj_e, local_nj_i, Nj_e, Nj_i ! Grids containing position in fourier space REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kxarray, kxarray_full REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kyarray, kyarray_full ! Kperp array depends on kx, ky, z (geometry), eo (even or odd zgrid) REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE, PUBLIC :: kparray REAL(dp), PUBLIC, PROTECTED :: deltakx, deltaky, kx_max, ky_max, kx_min, ky_min!, kp_max REAL(dp), PUBLIC, PROTECTED :: local_kxmax, local_kymax INTEGER, PUBLIC, PROTECTED :: ikxs, ikxe, ikys, ikye!, ikps, ikpe INTEGER, PUBLIC, PROTECTED :: ikx_0, iky_0, ikx_max, iky_max ! Indices of k-grid origin and max INTEGER, PUBLIC :: ikx, iky, ip, ij, ikp, pp2, eo ! counters LOGICAL, PUBLIC, PROTECTED :: contains_kx0 = .false. ! flag if the proc contains kx=0 index LOGICAL, PUBLIC, PROTECTED :: contains_ky0 = .false. ! flag if the proc contains ky=0 index LOGICAL, PUBLIC, PROTECTED :: contains_kymax = .false. ! flag if the proc contains kx=kxmax index LOGICAL, PUBLIC, PROTECTED :: contains_zmax = .false. ! flag if the proc contains z=pi-dz index LOGICAL, PUBLIC, PROTECTED :: contains_zmin = .false. ! flag if the proc contains z=-pi index LOGICAL, PUBLIC, PROTECTED :: SINGLE_KY = .false. ! to check if it is a single non 0 ky simulation ! Grid containing the polynomials degrees INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_e, parray_e_full INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_i, parray_i_full INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_e, jarray_e_full INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_i, jarray_i_full INTEGER, PUBLIC, PROTECTED :: ips_e,ipe_e, ijs_e,ije_e ! Start and end indices for pol. deg. INTEGER, PUBLIC, PROTECTED :: ips_i,ipe_i, ijs_i,ije_i INTEGER, PUBLIC, PROTECTED :: ipgs_e,ipge_e, ijgs_e,ijge_e ! Ghosts start and end indices INTEGER, PUBLIC, PROTECTED :: ipgs_i,ipge_i, ijgs_i,ijge_i INTEGER, PUBLIC, PROTECTED :: deltape, ip0_e, ip1_e, ip2_e, ip3_e ! Pgrid spacing and moment 0,1,2 index INTEGER, PUBLIC, PROTECTED :: deltapi, ip0_i, ip1_i, ip2_i, ip3_i LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip0_e, CONTAINS_ip0_i LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip1_e, CONTAINS_ip1_i LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip2_e, CONTAINS_ip2_i LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip3_e, CONTAINS_ip3_i LOGICAL, PUBLIC, PROTECTED :: SOLVE_POISSON, SOLVE_AMPERE INTEGER, PUBLIC, PROTECTED :: ij0_i, ij0_e ! Usefull inverse numbers REAL(dp), PUBLIC, PROTECTED :: inv_Nx, inv_Ny, inv_Nz ! Public Functions PUBLIC :: init_1Dgrid_distr PUBLIC :: set_pgrid, set_jgrid PUBLIC :: set_kxgrid, set_kygrid, set_zgrid PUBLIC :: grid_readinputs, grid_outputinputs PUBLIC :: bare, bari ! Precomputations real(dp), PUBLIC, PROTECTED :: pmaxe_dp, pmaxi_dp, jmaxe_dp,jmaxi_dp CONTAINS SUBROUTINE grid_readinputs ! Read the input parameters USE prec_const IMPLICIT NONE INTEGER :: lu_in = 90 ! File duplicated from STDIN NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, & - Nx, Lx, Nexc, Ny, Ly, Nz, Npol, SG + Nx, Lx, Ny, Ly, Nz, Npol, SG, Nexc READ(lu_in,grid) IF(Nz .EQ. 1) & ! overwrite SG option if Nz = 1 for safety of use SG = .FALSE. !! Compute the maximal degree of full GF moments set ! i.e. : all moments N_a^pj s.t. p+2j<=d are simulated (see GF closure) dmaxe = min(pmaxe,2*jmaxe+1) dmaxi = min(pmaxi,2*jmaxi+1) ! Usefull precomputations inv_Nx = 1._dp/REAL(Nx,dp) inv_Ny = 1._dp/REAL(Ny,dp) END SUBROUTINE grid_readinputs SUBROUTINE init_1Dgrid_distr ! write(*,*) Nx local_nky = (Ny/2+1)/num_procs_ky ! write(*,*) local_nkx local_nky_offset = rank_ky*local_nky if (rank_ky .EQ. num_procs_ky-1) local_nky = (Ny/2+1)-local_nky_offset END SUBROUTINE init_1Dgrid_distr SUBROUTINE set_pgrid USE prec_const USE model, ONLY: beta ! To know if we solve ampere or not and put odd p moments IMPLICIT NONE INTEGER :: ip, istart, iend, in ! If no parallel dim (Nz=1) and no EM effects (beta=0), the moment hierarchy !! is separable between odds and even P and since the energy is injected in !! P=0 and P=2 for density/temperature gradients there is no need of !! simulating the odd p which will only be damped. !! We define in this case a grid Parray = 0,2,4,...,Pmax i.e. deltap = 2 !! instead of 1 to spare computation IF((Nz .EQ. 1) .AND. (beta .EQ. 0._dp)) THEN deltape = 2; deltapi = 2; pp2 = 1; ! index p+2 is ip+1 ELSE deltape = 1; deltapi = 1; pp2 = 2; ! index p+2 is ip+2 ENDIF ! Total number of Hermite polynomials we will evolve total_np_e = (Pmaxe/deltape) + 1 total_np_i = (Pmaxi/deltapi) + 1 Np_e = total_np_e ! Reduced names (redundant) Np_i = total_np_i ! Build the full grids on process 0 to diagnose it without comm ALLOCATE(parray_e_full(1:total_np_e)) ALLOCATE(parray_i_full(1:total_np_i)) ! P DO ip = 1,total_np_e; parray_e_full(ip) = (ip-1)*deltape; END DO DO ip = 1,total_np_i; parray_i_full(ip) = (ip-1)*deltapi; END DO !! Parallel data distribution ! Local data distribution CALL decomp1D(total_np_e, num_procs_p, rank_p, ips_e, ipe_e) CALL decomp1D(total_np_i, num_procs_p, rank_p, ips_i, ipe_i) local_np_e = ipe_e - ips_e + 1 local_np_i = ipe_i - ips_i + 1 ! Ghosts boundaries ipgs_e = ips_e - 2/deltape; ipge_e = ipe_e + 2/deltape; ipgs_i = ips_i - 2/deltapi; ipge_i = ipe_i + 2/deltapi; ! List of shift and local numbers between the different processes (used in scatterv and gatherv) ALLOCATE(rcv_p_e (1:num_procs_p)) ALLOCATE(rcv_p_i (1:num_procs_p)) ALLOCATE(dsp_p_e (1:num_procs_p)) ALLOCATE(dsp_p_i (1:num_procs_p)) DO in = 0,num_procs_p-1 CALL decomp1D(total_np_e, num_procs_p, in, istart, iend) rcv_p_e(in+1) = iend-istart+1 dsp_p_e(in+1) = istart-1 CALL decomp1D(total_np_i, num_procs_p, in, istart, iend) rcv_p_i(in+1) = iend-istart+1 dsp_p_i(in+1) = istart-1 ENDDO !! local grid computations ! Flag to avoid unnecessary logical operations CONTAINS_ip0_e = .FALSE.; CONTAINS_ip1_e = .FALSE. CONTAINS_ip2_e = .FALSE.; CONTAINS_ip3_e = .FALSE. CONTAINS_ip0_i = .FALSE.; CONTAINS_ip1_i = .FALSE. CONTAINS_ip2_i = .FALSE.; CONTAINS_ip3_i = .FALSE. SOLVE_POISSON = .FALSE.; SOLVE_AMPERE = .FALSE. ALLOCATE(parray_e(ipgs_e:ipge_e)) ALLOCATE(parray_i(ipgs_i:ipge_i)) DO ip = ipgs_e,ipge_e parray_e(ip) = (ip-1)*deltape ! Storing indices of particular degrees for fluid moments computations SELECT CASE (parray_e(ip)) CASE(0); ip0_e = ip; CONTAINS_ip0_e = .TRUE. CASE(1); ip1_e = ip; CONTAINS_ip1_e = .TRUE. CASE(2); ip2_e = ip; CONTAINS_ip2_e = .TRUE. CASE(3); ip3_e = ip; CONTAINS_ip3_e = .TRUE. END SELECT END DO DO ip = ipgs_i,ipge_i parray_i(ip) = (ip-1)*deltapi ! Storing indices of particular degrees for fluid moments computations SELECT CASE (parray_i(ip)) CASE(0); ip0_i = ip; CONTAINS_ip0_i = .TRUE. CASE(1); ip1_i = ip; CONTAINS_ip1_i = .TRUE. CASE(2); ip2_i = ip; CONTAINS_ip2_i = .TRUE. CASE(3); ip3_i = ip; CONTAINS_ip3_i = .TRUE. END SELECT END DO IF(CONTAINS_ip0_e .AND. CONTAINS_ip0_i) SOLVE_POISSON = .TRUE. IF(CONTAINS_ip1_e .AND. CONTAINS_ip1_i) SOLVE_AMPERE = .TRUE. !DGGK operator uses moments at index p=2 (ip=3) for the p=0 term so the ! process that contains ip=1 MUST contain ip=3 as well for both e and i. IF(((ips_e .EQ. ip0_e) .OR. (ips_i .EQ. ip0_e)) .AND. ((ipe_e .LT. ip2_e) .OR. (ipe_i .LT. ip2_i)))& WRITE(*,*) "Warning : distribution along p may not work with DGGK" ! Precomputations pmaxe_dp = real(pmaxe,dp) pmaxi_dp = real(pmaxi,dp) ! Overwrite SOLVE_AMPERE flag if beta is zero IF(beta .EQ. 0._dp) THEN SOLVE_AMPERE = .FALSE. ENDIF END SUBROUTINE set_pgrid SUBROUTINE set_jgrid USE prec_const IMPLICIT NONE INTEGER :: ij ! Total number of J degrees Nj_e = jmaxe+1 Nj_i = jmaxi+1 ! Build the full grids on process 0 to diagnose it without comm ALLOCATE(jarray_e_full(1:jmaxe+1)) ALLOCATE(jarray_i_full(1:jmaxi+1)) ! J DO ij = 1,jmaxe+1; jarray_e_full(ij) = (ij-1); END DO DO ij = 1,jmaxi+1; jarray_i_full(ij) = (ij-1); END DO ! Local data ijs_e = 1; ije_e = jmaxe + 1 ijs_i = 1; ije_i = jmaxi + 1 ! Ghosts boundaries ijgs_e = ijs_e - 1; ijge_e = ije_e + 1; ijgs_i = ijs_i - 1; ijge_i = ije_i + 1; ! Local number of J local_nj_e = ijge_e - ijgs_e + 1 local_nj_i = ijge_i - ijgs_i + 1 ALLOCATE(jarray_e(ijgs_e:ijge_e)) ALLOCATE(jarray_i(ijgs_i:ijge_i)) DO ij = ijgs_e,ijge_e; jarray_e(ij) = ij-1; END DO DO ij = ijgs_i,ijge_i; jarray_i(ij) = ij-1; END DO ! Precomputations maxj = MAX(jmaxi, jmaxe) jmaxe_dp = real(jmaxe,dp) jmaxi_dp = real(jmaxi,dp) ! j=0 indices DO ij = ijs_e,ije_e; IF(jarray_e(ij) .EQ. 0) ij0_e = ij; END DO DO ij = ijs_i,ije_i; IF(jarray_i(ij) .EQ. 0) ij0_i = ij; END DO END SUBROUTINE set_jgrid SUBROUTINE set_kygrid USE prec_const USE model, ONLY: LINEARITY, N_HD IMPLICIT NONE INTEGER :: in, istart, iend Nky = Ny/2+1 ! Defined only on positive kx since fields are real ! Grid spacings IF (Ny .EQ. 1) THEN deltaky = 2._dp*PI/Ly ky_max = deltaky ky_min = deltaky ELSE deltaky = 2._dp*PI/Ly ky_max = (Nky-1)*deltaky ky_min = deltaky ENDIF ! Build the full grids on process 0 to diagnose it without comm ALLOCATE(kyarray_full(1:Nky)) DO iky = 1,Nky kyarray_full(iky) = REAL(iky-1,dp) * deltaky END DO !! Parallel distribution ikys = local_nky_offset + 1 ikye = ikys + local_nky - 1 ALLOCATE(kyarray(ikys:ikye)) local_kymax = 0._dp ! List of shift and local numbers between the different processes (used in scatterv and gatherv) ALLOCATE(counts_nky (1:num_procs_ky)) ALLOCATE(displs_nky (1:num_procs_ky)) DO in = 0,num_procs_ky-1 CALL decomp1D(Nky, num_procs_ky, in, istart, iend) counts_nky(in+1) = iend-istart+1 displs_nky(in+1) = istart-1 ENDDO ! Creating a grid ordered as dk*(0 1 2 3) DO iky = ikys,ikye IF(Ny .EQ. 1) THEN kyarray(iky) = deltaky kyarray_full(iky) = deltaky SINGLE_KY = .TRUE. ELSE kyarray(iky) = REAL(iky-1,dp) * deltaky ENDIF ! Finding kx=0 IF (kyarray(iky) .EQ. 0) THEN iky_0 = iky contains_ky0 = .true. ENDIF ! Finding local kxmax value IF (ABS(kyarray(iky)) .GT. local_kymax) THEN local_kymax = ABS(kyarray(iky)) ENDIF ! Finding kxmax idx IF (kyarray(iky) .EQ. ky_max) THEN iky_max = iky contains_kymax = .true. ENDIF END DO ! Orszag 2/3 filter two_third_kymax = 2._dp/3._dp*deltaky*(Nky-1) ! For hyperdiffusion IF(LINEARITY.EQ.'linear') THEN diff_ky_coeff= (1._dp/ky_max)**N_HD ELSE diff_ky_coeff= (1._dp/two_third_kymax)**N_HD ENDIF ALLOCATE(AA_y(ikys:ikye)) DO iky = ikys,ikye IF ( (kyarray(iky) .LT. two_third_kymax) .OR. (LINEARITY .EQ. 'linear')) THEN AA_y(iky) = 1._dp; ELSE AA_y(iky) = 0._dp; ENDIF END DO END SUBROUTINE set_kygrid SUBROUTINE set_kxgrid(shear) USE prec_const USE model, ONLY: LINEARITY, N_HD IMPLICIT NONE REAL(dp), INTENT(IN) :: shear REAL :: Lx_adapted IF(shear .GT. 0) THEN IF(my_id.EQ.0) write(*,*) 'Magnetic shear detected: set up sheared kx grid..' ! mininal size of box in x to respect dkx = 2pi shear dky Lx_adapted = Ly/(2._dp*pi*shear*Npol) ! Put Nexc to 0 so that it is computed from a target value Lx IF(Nexc .EQ. 0) THEN Nexc = CEILING(0.9 * Lx/Lx_adapted) IF(my_id.EQ.0) write(*,*) 'Adapted Nexc =', Nexc ENDIF ! x length is adapted Lx = Lx_adapted*Nexc ENDIF Nkx = Nx; ! Local data ! Start and END indices of grid ikxs = 1 ikxe = Nkx local_nkx = ikxe - ikxs + 1 ALLOCATE(kxarray(ikxs:ikxe)) ALLOCATE(kxarray_full(1:Nkx)) IF (Nx .EQ. 1) THEN deltakx = 1._dp kxarray(1) = 0._dp ikx_0 = 1 contains_kx0 = .true. kx_max = 0._dp ikx_max = 1 kx_min = 0._dp kxarray_full(1) = 0._dp local_kxmax = 0._dp ELSE ! Build apprpopriate grid deltakx = 2._dp*PI/Lx IF(MODULO(Nkx,2) .EQ. 0) THEN ! Even number of Nkx (-2 -1 0 1 2 3) kx_max = (Nkx/2)*deltakx kx_min = -kx_max+deltakx ! Creating a grid ordered as dk*(0 1 2 3 -2 -1) local_kxmax = 0._dp DO ikx = ikxs,ikxe kxarray(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1)/real(Nkx))) if (ikx .EQ. Nx/2+1) kxarray(ikx) = -kxarray(ikx) ! Finding kx=0 IF (kxarray(ikx) .EQ. 0) THEN ikx_0 = ikx contains_kx0 = .true. ENDIF ! Finding local kxmax IF (ABS(kxarray(ikx)) .GT. local_kxmax) THEN local_kxmax = ABS(kxarray(ikx)) ENDIF ! Finding kxmax IF (kxarray(ikx) .EQ. kx_max) ikx_max = ikx END DO ! Build the full grids on process 0 to diagnose it without comm ! kx DO ikx = 1,Nkx kxarray_full(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1)/real(Nkx))) IF (ikx .EQ. Nx/2+1) kxarray_full(ikx) = -kxarray_full(ikx) END DO ELSE ! Odd number of kx (-2 -1 0 1 2) kx_max = (Nkx-1)/2*deltakx kx_min = -kx_max ! Creating a grid ordered as dk*(0 1 2 -2 -1) local_kxmax = 0._dp DO ikx = ikxs,ikxe IF(ikx .LE. (Nkx-1)/2+1) THEN kxarray(ikx) = deltakx*(ikx-1) ELSE kxarray(ikx) = deltakx*(ikx-Nkx-1) ENDIF ! Finding kx=0 IF (kxarray(ikx) .EQ. 0) THEN ikx_0 = ikx contains_kx0 = .true. ENDIF ! Finding local kxmax IF (ABS(kxarray(ikx)) .GT. local_kxmax) THEN local_kxmax = ABS(kxarray(ikx)) ENDIF ! Finding kxmax IF (kxarray(ikx) .EQ. kx_max) ikx_max = ikx END DO ! Build the full grids on process 0 to diagnose it without comm ! kx DO ikx = 1,Nkx IF(ikx .LE. (Nkx-1)/2+1) THEN kxarray_full(ikx) = deltakx*(ikx-1) ELSE kxarray_full(ikx) = deltakx*(ikx-Nkx-1) ENDIF END DO ENDIF ENDIF ! Orszag 2/3 filter two_third_kxmax = 2._dp/3._dp*kx_max; ! For hyperdiffusion IF(LINEARITY.EQ.'linear') THEN diff_kx_coeff= (1._dp/kx_max)**N_HD ELSE diff_kx_coeff= (1._dp/two_third_kxmax)**N_HD ENDIF ! Antialiasing filter ALLOCATE(AA_x(ikxs:ikxe)) DO ikx = ikxs,ikxe IF ( ((kxarray(ikx) .GT. -two_third_kxmax) .AND. & (kxarray(ikx) .LT. two_third_kxmax)) .OR. (LINEARITY .EQ. 'linear')) THEN AA_x(ikx) = 1._dp; ELSE AA_x(ikx) = 0._dp; ENDIF END DO END SUBROUTINE set_kxgrid SUBROUTINE set_zgrid USE prec_const USE model, ONLY: mu_z IMPLICIT NONE REAL :: grid_shift, Lz, zmax, zmin INTEGER :: istart, iend, in total_nz = Nz ! Length of the flux tube (in ballooning angle) Lz = 2_dp*pi*Npol ! Z stepping (#interval = #points since periodic) deltaz = Lz/REAL(Nz,dp) inv_deltaz = 1._dp/deltaz ! Parallel hyperdiffusion coefficient IF(mu_z .GT. 0) THEN diff_dz_coeff = REAL((imagu*deltaz/2._dp)**4) ! adaptive fourth derivative (~GENE) ELSE diff_dz_coeff = 1._dp ! non adaptive (positive sign to compensate mu_z neg) ENDIF IF (SG) THEN grid_shift = deltaz/2._dp ELSE grid_shift = 0._dp ENDIF ! Build the full grids on process 0 to diagnose it without comm ALLOCATE(zarray_full(1:Nz)) IF (Nz .EQ. 1) Npol = 0 zmax = 0; zmin = 0; DO iz = 1,total_nz ! z in [-pi pi-dz] x Npol zarray_full(iz) = REAL(iz-1,dp)*deltaz - Lz/2._dp IF(zarray_full(iz) .GT. zmax) zmax = zarray_full(iz) IF(zarray_full(iz) .LT. zmin) zmin = zarray_full(iz) END DO !! Parallel data distribution - IF( num_procs_z .GT. 1) stop '>>STOPPED<< Cannot have multiple core in z-direction (Nz = 1)' + IF( (Nz .EQ. 1) .AND. (num_procs_z .GT. 1) ) & + stop '>>STOPPED<< Cannot have multiple core in z-direction (Nz = 1)' ! Local data distribution CALL decomp1D(total_nz, num_procs_z, rank_z, izs, ize) local_nz = ize - izs + 1 ! Ghosts boundaries (depend on the order of z operators) IF(Nz .EQ. 1) THEN izgs = izs; izge = ize; zarray_full(izs) = 0; ELSEIF(Nz .GE. 4) THEN izgs = izs - 2; izge = ize + 2; ELSE ERROR STOP 'Error stop: Nz is not appropriate!!' ENDIF ! List of shift and local numbers between the different processes (used in scatterv and gatherv) ALLOCATE(counts_nz (1:num_procs_z)) ALLOCATE(displs_nz (1:num_procs_z)) DO in = 0,num_procs_z-1 CALL decomp1D(total_nz, num_procs_z, in, istart, iend) counts_nz(in+1) = iend-istart+1 displs_nz(in+1) = istart-1 ENDDO ! Local z array ALLOCATE(zarray(izgs:izge,0:1)) DO iz = izgs,izge IF(iz .EQ. 0) THEN zarray(iz,0) = zarray_full(total_nz) zarray(iz,1) = zarray_full(total_nz) + grid_shift ELSEIF(iz .EQ. -1) THEN zarray(iz,0) = zarray_full(total_nz-1) zarray(iz,1) = zarray_full(total_nz-1) + grid_shift ELSEIF(iz .EQ. total_nz + 1) THEN zarray(iz,0) = zarray_full(1) zarray(iz,1) = zarray_full(1) + grid_shift ELSEIF(iz .EQ. total_nz + 2) THEN zarray(iz,0) = zarray_full(2) zarray(iz,1) = zarray_full(2) + grid_shift ELSE zarray(iz,0) = zarray_full(iz) zarray(iz,1) = zarray_full(iz) + grid_shift ENDIF ENDDO IF(abs(zarray(izs,0) - zmin) .LT. EPSILON(zmin)) & contains_zmin = .TRUE. IF(abs(zarray(ize,0) - zmax) .LT. EPSILON(zmax)) & contains_zmax = .TRUE. ! Weitghs for Simpson rule ALLOCATE(zweights_SR(izs:ize)) DO iz = izs,ize IF(MODULO(iz,2) .EQ. 1) THEN ! odd iz zweights_SR(iz) = 2._dp ELSE ! even iz zweights_SR(iz) = 4._dp ENDIF ENDDO END SUBROUTINE set_zgrid SUBROUTINE grid_outputinputs(fidres, str) ! Write the input parameters to the results_xx.h5 file USE futils, ONLY: attach USE prec_const IMPLICIT NONE INTEGER, INTENT(in) :: fidres CHARACTER(len=256), INTENT(in) :: str CALL attach(fidres, TRIM(str), "pmaxe", pmaxe) CALL attach(fidres, TRIM(str), "jmaxe", jmaxe) CALL attach(fidres, TRIM(str), "pmaxi", pmaxi) CALL attach(fidres, TRIM(str), "jmaxi", jmaxi) CALL attach(fidres, TRIM(str), "Nx", Nx) CALL attach(fidres, TRIM(str), "Lx", Lx) CALL attach(fidres, TRIM(str), "Nexc", Nexc) CALL attach(fidres, TRIM(str), "Ny", Ny) CALL attach(fidres, TRIM(str), "Ly", Ly) CALL attach(fidres, TRIM(str), "Nz", Nz) CALL attach(fidres, TRIM(str), "Nkx", Nkx) CALL attach(fidres, TRIM(str), "Lkx", Lkx) CALL attach(fidres, TRIM(str), "Nky", Nky) CALL attach(fidres, TRIM(str), "Lky", Lky) CALL attach(fidres, TRIM(str), "SG", SG) END SUBROUTINE grid_outputinputs FUNCTION bare(p_,j_) IMPLICIT NONE INTEGER :: bare, p_, j_ bare = (jmaxe+1)*p_ + j_ + 1 END FUNCTION FUNCTION bari(p_,j_) IMPLICIT NONE INTEGER :: bari, p_, j_ bari = (jmaxi+1)*p_ + j_ + 1 END FUNCTION SUBROUTINE decomp1D( n, numprocs, myid, s, e ) INTEGER :: n, numprocs, myid, s, e INTEGER :: nlocal INTEGER :: deficit nlocal = n / numprocs s = myid * nlocal + 1 deficit = MOD(n,numprocs) s = s + MIN(myid,deficit) IF (myid .LT. deficit) nlocal = nlocal + 1 e = s + nlocal - 1 IF (e .GT. n .OR. myid .EQ. numprocs-1) e = n END SUBROUTINE decomp1D END MODULE grid diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90 index 25e7f57..d866c94 100644 --- a/src/moments_eq_rhs_mod.F90 +++ b/src/moments_eq_rhs_mod.F90 @@ -1,409 +1,409 @@ MODULE moments_eq_rhs IMPLICIT NONE PUBLIC :: compute_moments_eq_rhs CONTAINS SUBROUTINE compute_moments_eq_rhs USE model, only: KIN_E IMPLICIT NONE IF(KIN_E) CALL moments_eq_rhs_e CALL moments_eq_rhs_i !! BETA TESTING !! IF(.false.) CALL add_Maxwellian_background_terms END SUBROUTINE compute_moments_eq_rhs !_____________________________________________________________________________! !_____________________________________________________________________________! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! Electrons moments RHS !!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !_____________________________________________________________________________! SUBROUTINE moments_eq_rhs_e USE basic USE time_integration USE array, ONLY: xnepj, xnepp2j, xnepm2j, xnepjp1, xnepjm1, xnepp1j, xnepm1j,& ynepp1j, ynepp1jm1, ynepm1j, ynepm1jm1, & znepm1j, znepm1jp1, znepm1jm1, & xphij_e, xphijp1_e, xphijm1_e, xpsij_e, xpsijp1_e, xpsijm1_e,& nadiab_moments_e, ddz_nepj, interp_nepj, moments_rhs_e, Sepj,& TColl_e, ddzND_nepj, kernel_e USE fields, ONLY: phi, psi, moments_e USE grid USE model USE prec_const USE collision USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatB USE calculus, ONLY : interp_z, grad_z, grad_z2 IMPLICIT NONE INTEGER :: p_int, j_int ! loops indices and polynom. degrees REAL(dp) :: kx, ky, kperp2 COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1 ! Terms from b x gradB and drives COMPLEX(dp) :: Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi, Tpsi COMPLEX(dp) :: Unepm1j, Unepm1jp1, Unepm1jm1 ! Terms from mirror force with adiab moments COMPLEX(dp) :: i_kx,i_ky,phikykxz, psikykxz ! Measuring execution time CALL cpu_time(t0_rhs) ! Spatial loops zloope : DO iz = izs,ize kxloope : DO ikx = ikxs,ikxe kx = kxarray(ikx) ! radial wavevector i_kx = imagu * kx ! toroidal derivative kyloope : DO iky = ikys,ikye ky = kyarray(iky) ! toroidal wavevector i_ky = imagu * ky ! toroidal derivative phikykxz = phi(iky,ikx,iz)! tmp phi value psikykxz = psi(iky,ikx,iz)! tmp psi value ! Kinetic loops jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxi+1 j_int = jarray_e(ij) ploope : DO ip = ips_e, ipe_e ! Hermite loop p_int = parray_e(ip) ! Hermite degree eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid kperp2= kparray(iky,ikx,iz,eo)**2 IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN !! Compute moments mixing terms Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp ! Perpendicular dynamic ! term propto n_e^{p,j} Tnepj = xnepj(ip,ij)* nadiab_moments_e(ip,ij,iky,ikx,iz) ! term propto n_e^{p+2,j} Tnepp2j = xnepp2j(ip) * nadiab_moments_e(ip+pp2,ij,iky,ikx,iz) ! term propto n_e^{p-2,j} Tnepm2j = xnepm2j(ip) * nadiab_moments_e(ip-pp2,ij,iky,ikx,iz) ! term propto n_e^{p,j+1} Tnepjp1 = xnepjp1(ij) * nadiab_moments_e(ip,ij+1,iky,ikx,iz) ! term propto n_e^{p,j-1} Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,iky,ikx,iz) ! Tperp Tperp = Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1 ! Parallel dynamic ! ddz derivative for Landau damping term Tpar = xnepp1j(ip) * ddz_nepj(ip+1,ij,iky,ikx,iz) & + xnepm1j(ip) * ddz_nepj(ip-1,ij,iky,ikx,iz) ! Mirror terms Tnepp1j = ynepp1j (ip,ij) * interp_nepj(ip+1,ij ,iky,ikx,iz) Tnepp1jm1 = ynepp1jm1(ip,ij) * interp_nepj(ip+1,ij-1,iky,ikx,iz) Tnepm1j = ynepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz) Tnepm1jm1 = ynepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz) ! Trapping terms Unepm1j = znepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz) Unepm1jp1 = znepm1jp1(ip,ij) * interp_nepj(ip-1,ij+1,iky,ikx,iz) Unepm1jm1 = znepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz) Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + Unepm1j + Unepm1jp1 + Unepm1jm1 !! Electrical potential term IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 Tphi = (xphij_e (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) & + xphijp1_e(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) & + xphijm1_e(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*phikykxz ELSE Tphi = 0._dp ENDIF !! Vector potential term IF ( (p_int .LE. 3) .AND. (p_int .GE. 1) ) THEN ! Kronecker p1 or p3 Tpsi = (xpsij_e (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) & + xpsijp1_e(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) & + xpsijm1_e(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*psikykxz ELSE Tpsi = 0._dp ENDIF !! Sum of all RHS terms moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = & ! Perpendicular magnetic gradient/curvature effects -imagu*Ckxky(iky,ikx,iz,eo)*hatB(iz,eo) * Tperp& ! Perpendicular pressure effects -i_ky*beta*dpdx * Tperp& ! Parallel coupling (Landau Damping) -gradz_coeff(iz,eo) * Tpar & ! Mirror term (parallel magnetic gradient) -gradzB(iz,eo)*gradz_coeff(iz,eo) * Tmir& ! Drives (density + temperature gradients) -i_ky * (Tphi - Tpsi) & ! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose) -mu_x*diff_kx_coeff*kx**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) & -mu_y*diff_ky_coeff*ky**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) & ! Numerical parallel hyperdiffusion "mu_z*ddz**4" see Pueschel 2010 (eq 25) -mu_z*diff_dz_coeff*ddzND_nepj(ip,ij,iky,ikx,iz) & ! Collision term +TColl_e(ip,ij,iky,ikx,iz) & ! Nonlinear term -hatB_NL(iz,eo) * Sepj(ip,ij,iky,ikx,iz) - IF( (ip-4 .GT. 0) .AND. (num_procs_p .EQ. 1) ) & - ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33) - ! (not used often so not parallelized) - moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = & - moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) & - + mu_p * moments_e(ip-4,ij,iky,ikx,iz,updatetlevel) + ! IF( (ip-4 .GT. 0) .AND. (num_procs_p .EQ. 1) ) & + ! ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33) + ! ! (not used often so not parallelized) + ! moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = & + ! moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) & + ! + mu_p * moments_e(ip-4,ij,iky,ikx,iz,updatetlevel) ELSE moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp ENDIF END DO ploope END DO jloope END DO kyloope END DO kxloope END DO zloope ! Execution time end CALL cpu_time(t1_rhs) tc_rhs = tc_rhs + (t1_rhs-t0_rhs) END SUBROUTINE moments_eq_rhs_e !_____________________________________________________________________________! !_____________________________________________________________________________! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!! Ions moments RHS !!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !_____________________________________________________________________________! SUBROUTINE moments_eq_rhs_i USE basic USE time_integration, ONLY: updatetlevel USE array, ONLY: xnipj, xnipp2j, xnipm2j, xnipjp1, xnipjm1, xnipp1j, xnipm1j,& ynipp1j, ynipp1jm1, ynipm1j, ynipm1jm1, & znipm1j, znipm1jp1, znipm1jm1, & xphij_i, xphijp1_i, xphijm1_i, xpsij_i, xpsijp1_i, xpsijm1_i,& nadiab_moments_i, ddz_nipj, interp_nipj, moments_rhs_i, Sipj,& TColl_i, ddzND_nipj, kernel_i USE fields, ONLY: phi, psi, moments_i USE grid USE model USE prec_const USE collision USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatB USE calculus, ONLY : interp_z, grad_z, grad_z2 IMPLICIT NONE INTEGER :: p_int, j_int ! loops indices and polynom. degrees REAL(dp) :: kx, ky, kperp2 COMPLEX(dp) :: Tnipj, Tnipp2j, Tnipm2j, Tnipjp1, Tnipjm1 COMPLEX(dp) :: Tnipp1j, Tnipm1j, Tnipp1jm1, Tnipm1jm1 ! Terms from mirror force with non adiab moments COMPLEX(dp) :: Unipm1j, Unipm1jp1, Unipm1jm1 ! Terms from mirror force with adiab moments COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi, Tpsi COMPLEX(dp) :: i_kx, i_ky, phikykxz, psikykxz ! Measuring execution time CALL cpu_time(t0_rhs) ! Spatial loops zloopi : DO iz = izs,ize kxloopi : DO ikx = ikxs,ikxe kx = kxarray(ikx) ! radial wavevector i_kx = imagu * kx ! toroidal derivative kyloopi : DO iky = ikys,ikye ky = kyarray(iky) ! toroidal wavevector i_ky = imagu * ky ! toroidal derivative phikykxz = phi(iky,ikx,iz)! tmp phi value psikykxz = psi(iky,ikx,iz)! tmp phi value ! Kinetic loops jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1 j_int = jarray_i(ij) ploopi : DO ip = ips_i, ipe_i ! Hermite loop p_int = parray_i(ip) ! Hermite degree eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid kperp2= kparray(iky,ikx,iz,eo)**2 IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN !! Compute moments mixing terms Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp ! Perpendicular dynamic ! term propto n_i^{p,j} Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip ,ij ,iky,ikx,iz) ! term propto n_i^{p+2,j} Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij ,iky,ikx,iz) ! term propto n_i^{p-2,j} Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij ,iky,ikx,iz) ! term propto n_i^{p,j+1} Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip ,ij+1,iky,ikx,iz) ! term propto n_i^{p,j-1} Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip ,ij-1,iky,ikx,iz) ! Tperp Tperp = Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1 ! Parallel dynamic ! ddz derivative for Landau damping term Tpar = xnipp1j(ip) * ddz_nipj(ip+1,ij,iky,ikx,iz) & + xnipm1j(ip) * ddz_nipj(ip-1,ij,iky,ikx,iz) ! Mirror terms Tnipp1j = ynipp1j (ip,ij) * interp_nipj(ip+1,ij ,iky,ikx,iz) Tnipp1jm1 = ynipp1jm1(ip,ij) * interp_nipj(ip+1,ij-1,iky,ikx,iz) Tnipm1j = ynipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz) Tnipm1jm1 = ynipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz) ! Trapping terms Unipm1j = znipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz) Unipm1jp1 = znipm1jp1(ip,ij) * interp_nipj(ip-1,ij+1,iky,ikx,iz) Unipm1jm1 = znipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz) Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + Unipm1j + Unipm1jp1 + Unipm1jm1 !! Electrical potential term IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 Tphi = (xphij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) & + xphijp1_i(ip,ij)*kernel_i(ij+1,iky,ikx,iz,eo) & + xphijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*phikykxz ELSE Tphi = 0._dp ENDIF !! Vector potential term IF ( (p_int .LE. 3) .AND. (p_int .GE. 1) ) THEN ! Kronecker p1 or p3 Tpsi = (xpsij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) & + xpsijp1_i(ip,ij)*kernel_i(ij+1,iky,ikx,iz,eo) & + xpsijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*psikykxz ELSE Tpsi = 0._dp ENDIF !! Sum of all RHS terms moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = & ! Perpendicular magnetic gradient/curvature effects -imagu*Ckxky(iky,ikx,iz,eo)*hatB(iz,eo) * Tperp & ! Perpendicular pressure effects -i_ky*beta*dpdx * Tperp& ! Parallel coupling (Landau damping) -gradz_coeff(iz,eo) * Tpar & ! Mirror term (parallel magnetic gradient) -gradzB(iz,eo)*gradz_coeff(iz,eo) * Tmir & ! Drives (density + temperature gradients) -i_ky * (Tphi - Tpsi) & ! Numerical hyperdiffusion (totally artificial, for stability purpose) -mu_x*diff_kx_coeff*kx**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) & -mu_y*diff_ky_coeff*ky**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) & ! Numerical parallel hyperdiffusion "mu_z*ddz**4" -mu_z*diff_dz_coeff*ddzND_nipj(ip,ij,iky,ikx,iz) & ! Collision term +TColl_i(ip,ij,iky,ikx,iz)& ! Nonlinear term with a (gxx*gxy - gxy**2)^1/2 factor -hatB_NL(iz,eo) * Sipj(ip,ij,iky,ikx,iz) - IF( (ip-4 .GT. 0) .AND. (num_procs_p .EQ. 1) ) & - ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33) - ! (not used often so not parallelized) - moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = & - moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) & - + mu_p * moments_i(ip-4,ij,iky,ikx,iz,updatetlevel) + ! IF( (ip-4 .GT. 0) .AND. (num_procs_p .EQ. 1) ) & + ! ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33) + ! ! (not used often so not parallelized) + ! moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = & + ! moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) & + ! + mu_p * moments_i(ip-4,ij,iky,ikx,iz,updatetlevel) ELSE moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp ENDIF END DO ploopi END DO jloopi END DO kyloopi END DO kxloopi END DO zloopi ! Execution time end CALL cpu_time(t1_rhs) tc_rhs = tc_rhs + (t1_rhs-t0_rhs) END SUBROUTINE moments_eq_rhs_i SUBROUTINE add_Maxwellian_background_terms ! This routine is meant to add the terms rising from the magnetic operator, ! i.e. (B x GradB) Grad, applied on the background Maxwellian distribution ! (x_a + spar^2)(b x GradB) GradFaM ! It gives birth to kx=ky=0 sources terms (averages) that hit moments 00, 20, ! 40, 01,02, 21 with background gradient dependences. USE prec_const USE time_integration, ONLY : updatetlevel USE model, ONLY: taue_qe, taui_qi, k_Ni, k_Ne, k_Ti, k_Te, KIN_E USE array, ONLY: moments_rhs_e, moments_rhs_i USE grid, ONLY: contains_kx0, contains_ky0, ikx_0, iky_0,& ips_e,ipe_e,ijs_e,ije_e,ips_i,ipe_i,ijs_i,ije_i,& zarray, izs,ize,& ip,ij IMPLICIT NONE real(dp), DIMENSION(izs:ize) :: sinz sinz(izs:ize) = SIN(zarray(izs:ize,0)) IF(contains_kx0 .AND. contains_ky0) THEN IF(KIN_E) THEN DO ip = ips_e,ipe_e DO ij = ijs_e,ije_e SELECT CASE(ij-1) CASE(0) ! j = 0 SELECT CASE (ip-1) CASE(0) ! Na00 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taue_qe * sinz(izs:ize) * (1.5_dp*k_Ne - 1.125_dp*k_Te) CASE(2) ! Na20 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*k_Ne - 2.75_dp*k_Te) CASE(4) ! Na40 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*k_Te END SELECT CASE(1) ! j = 1 SELECT CASE (ip-1) CASE(0) ! Na01 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& -taue_qe * sinz(izs:ize) * (k_Ne + 3.5_dp*k_Te) CASE(2) ! Na21 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& -taue_qe * sinz(izs:ize) * SQRT2*k_Te END SELECT CASE(2) ! j = 2 SELECT CASE (ip-1) CASE(0) ! Na02 term moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taue_qe * sinz(izs:ize) * 2._dp*k_Te END SELECT END SELECT ENDDO ENDDO ENDIF DO ip = ips_i,ipe_i DO ij = ijs_i,ije_i SELECT CASE(ij-1) CASE(0) ! j = 0 SELECT CASE (ip-1) CASE(0) ! Na00 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taui_qi * sinz(izs:ize) * (1.5_dp*k_Ni - 1.125_dp*k_Ti) CASE(2) ! Na20 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*k_Ni - 2.75_dp*k_Ti) CASE(4) ! Na40 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*k_Ti END SELECT CASE(1) ! j = 1 SELECT CASE (ip-1) CASE(0) ! Na01 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& -taui_qi * sinz(izs:ize) * (k_Ni + 3.5_dp*k_Ti) CASE(2) ! Na21 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& -taui_qi * sinz(izs:ize) * SQRT2*k_Ti END SELECT CASE(2) ! j = 2 SELECT CASE (ip-1) CASE(0) ! Na02 term moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& +taui_qi * sinz(izs:ize) * 2._dp*k_Ti END SELECT END SELECT ENDDO ENDDO ENDIF END SUBROUTINE END MODULE moments_eq_rhs diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90 index 7203bdf..12b8291 100644 --- a/src/time_integration_mod.F90 +++ b/src/time_integration_mod.F90 @@ -1,226 +1,294 @@ MODULE time_integration USE prec_const IMPLICIT NONE PRIVATE INTEGER, PUBLIC, PROTECTED :: ntimelevel=4 ! Total number of time levels required by the numerical scheme INTEGER, PUBLIC, PROTECTED :: updatetlevel ! Current time level to be updated real(dp),PUBLIC,PROTECTED,DIMENSION(:,:),ALLOCATABLE :: A_E,A_I real(dp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: b_E,b_Es,b_I real(dp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: c_E,c_I !Coeff for Expl/Implic time integration in case of time dependent RHS (i.e. dy/dt = f(y,t)) see Baptiste Frei CSE Rapport 06/17 character(len=10),PUBLIC,PROTECTED :: numerical_scheme='RK4' PUBLIC :: set_updatetlevel, time_integration_readinputs, time_integration_outputinputs CONTAINS SUBROUTINE set_updatetlevel(new_updatetlevel) INTEGER, INTENT(in) :: new_updatetlevel updatetlevel = new_updatetlevel END SUBROUTINE set_updatetlevel SUBROUTINE time_integration_readinputs ! Read the input parameters USE prec_const USE basic, ONLY : lu_in IMPLICIT NONE NAMELIST /TIME_INTEGRATION_PAR/ numerical_scheme READ(lu_in,time_integration_par) CALL set_numerical_scheme END SUBROUTINE time_integration_readinputs SUBROUTINE time_integration_outputinputs(fidres, str) ! Write the input parameters to the results_xx.h5 file USE prec_const USE futils, ONLY: attach IMPLICIT NONE INTEGER, INTENT(in) :: fidres CHARACTER(len=256), INTENT(in) :: str CALL attach(fidres, TRIM(str), "numerical_scheme", numerical_scheme) END SUBROUTINE time_integration_outputinputs SUBROUTINE set_numerical_scheme ! Initialize Butcher coefficient of numerical_scheme use basic IMPLICIT NONE SELECT CASE (numerical_scheme) CASE ('RK4') CALL RK4 + CASE ('SSPx_RK3') + CALL SSPx_RK3 + CASE ('SSPx_RK2') + CALL SSPx_RK2 CASE ('DOPRI5') - CALL DOPRI5 - CASE ('DOPRI5_ADAPT') - CALL DOPRI5_ADAPT + CALL DOPRI5 CASE DEFAULT IF (my_id .EQ. 0) WRITE(*,*) 'Cannot initialize time integration scheme. Name invalid.' END SELECT IF (my_id .EQ. 0) WRITE(*,*) " Time integration with ", numerical_scheme END SUBROUTINE set_numerical_scheme SUBROUTINE RK4 ! Butcher coeff for RK4 (default) USE basic USE prec_const IMPLICIT NONE INTEGER,PARAMETER :: nbstep = 4 CALL allocate_array_dp1(c_E,1,nbstep) CALL allocate_array_dp1(b_E,1,nbstep) CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep) ntimelevel = 4 c_E(1) = 0.0_dp c_E(2) = 1.0_dp/2.0_dp c_E(3) = 1.0_dp/2.0_dp c_E(4) = 1.0_dp b_E(1) = 1.0_dp/6.0_dp b_E(2) = 1.0_dp/3.0_dp b_E(3) = 1.0_dp/3.0_dp b_E(4) = 1.0_dp/6.0_dp A_E(2,1) = 1.0_dp/2.0_dp A_E(3,2) = 1.0_dp/2.0_dp A_E(4,3) = 1.0_dp END SUBROUTINE RK4 + SUBROUTINE SSPx_RK3 + ! Butcher coeff for modified strong stability preserving RK3 + ! used in GX (Hammett 2022, Mandell et al. 2022) + USE basic + USE prec_const + IMPLICIT NONE + INTEGER,PARAMETER :: nbstep = 3 + REAL(dp) :: a1, a2, a3, w1, w2, w3 + a1 = (1._dp/6._dp)**(1._dp/3._dp)! (1/6)^(1/3) + ! a1 = 0.5503212081491044571635029569733887910843_dp ! (1/6)^(1/3) + a2 = a1 + a3 = a1 + w1 = 0.5_dp*(-1._dp + SQRT( 9._dp - 2._dp * 6._dp**(2._dp/3._dp))) ! (-1 + sqrt(9-2*6^(2/3)))/2 + ! w1 = 0.2739744023885328783052273138309828937054_dp ! (sqrt(9-2*6^(2/3))-1)/2 + w2 = 0.5_dp*(-1._dp + 6._dp**(2._dp/3._dp) - SQRT(9._dp - 2._dp * 6._dp**(2._dp/3._dp))) ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2 + ! w2 = 0.3769892220587804931852815570891834795475_dp ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2 + w3 = 1._dp/a1 - w2 * (1._dp + w1) + ! w3 = 1.3368459739528868457369981115334667265415_dp + + CALL allocate_array_dp1(c_E,1,nbstep) + CALL allocate_array_dp1(b_E,1,nbstep) + CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep) + + ntimelevel = 3 + + c_E(1) = 0.0_dp + c_E(2) = 1.0_dp/2.0_dp + c_E(3) = 1.0_dp/2.0_dp + + b_E(1) = a1 * (w1*w2 + w3) + b_E(2) = a2 * w2 + b_E(3) = a3 + + A_E(2,1) = a1 + A_E(3,1) = a1 * w1 + A_E(3,2) = a2 + + END SUBROUTINE SSPx_RK3 + + SUBROUTINE SSPx_RK2 + ! Butcher coeff for modified strong stability preserving RK2 + ! used in GX (Hammett 2022, Mandell et al. 2022) + + USE basic + USE prec_const + IMPLICIT NONE + INTEGER,PARAMETER :: nbstep = 2 + REAL(dp) :: alpha, beta + alpha = 1._dp/SQRT(2._dp) + beta = SQRT(2._dp) - 1._dp + + CALL allocate_array_dp1(c_E,1,nbstep) + CALL allocate_array_dp1(b_E,1,nbstep) + CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep) + + ntimelevel = 2 + + c_E(1) = 0.0_dp + c_E(2) = 1.0_dp/2.0_dp + + b_E(1) = alpha*beta/2._dp + b_E(2) = alpha/2._dp + + A_E(2,1) = alpha + + END SUBROUTINE SSPx_RK2 SUBROUTINE DOPRI5 ! Butcher coeff for DOPRI5 --> Stiffness detection ! DOPRI5 used for stiffness detection. ! 5 order method/7 stages USE basic IMPLICIT NONE INTEGER,PARAMETER :: nbstep =7 CALL allocate_array_dp1(c_E,1,nbstep) CALL allocate_array_dp1(b_E,1,nbstep) CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep) ntimelevel = 7 ! c_E(1) = 0._dp c_E(2) = 1.0_dp/5.0_dp c_E(3) = 3.0_dp /10.0_dp c_E(4) = 4.0_dp/5.0_dp c_E(5) = 8.0_dp/9.0_dp c_E(6) = 1.0_dp c_E(7) = 1.0_dp ! A_E(2,1) = 1.0_dp/5.0_dp A_E(3,1) = 3.0_dp/40.0_dp A_E(3,2) = 9.0_dp/40.0_dp A_E(4,1) = 44.0_dp/45.0_dp A_E(4,2) = -56.0_dp/15.0_dp A_E(4,3) = 32.0_dp/9.0_dp A_E(5,1 ) = 19372.0_dp/6561.0_dp A_E(5,2) = -25360.0_dp/2187.0_dp A_E(5,3) = 64448.0_dp/6561.0_dp A_E(5,4) = -212.0_dp/729.0_dp A_E(6,1) = 9017.0_dp/3168.0_dp A_E(6,2)= -355.0_dp/33.0_dp A_E(6,3) = 46732.0_dp/5247.0_dp A_E(6,4) = 49.0_dp/176.0_dp A_E(6,5) = -5103.0_dp/18656.0_dp A_E(7,1) = 35.0_dp/384.0_dp A_E(7,3) = 500.0_dp/1113.0_dp A_E(7,4) = 125.0_dp/192.0_dp A_E(7,5) = -2187.0_dp/6784.0_dp A_E(7,6) = 11.0_dp/84.0_dp ! b_E(1) = 35.0_dp/384.0_dp b_E(2) = 0._dp b_E(3) = 500.0_dp/1113.0_dp b_E(4) = 125.0_dp/192.0_dp b_E(5) = -2187.0_dp/6784.0_dp b_E(6) = 11.0_dp/84.0_dp b_E(7) = 0._dp ! END SUBROUTINE DOPRI5 SUBROUTINE DOPRI5_ADAPT ! Butcher coeff for DOPRI5 --> Stiffness detection ! DOPRI5 used for stiffness detection. ! 5 order method/7 stages USE basic IMPLICIT NONE INTEGER,PARAMETER :: nbstep =7 CALL allocate_array_dp1(c_E,1,nbstep) CALL allocate_array_dp1(b_E,1,nbstep) CALL allocate_array_dp1(b_Es,1,nbstep) CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep) ntimelevel = 7 ! c_E(1) = 0._dp c_E(2) = 1.0_dp/5.0_dp c_E(3) = 3.0_dp /10.0_dp c_E(4) = 4.0_dp/5.0_dp c_E(5) = 8.0_dp/9.0_dp c_E(6) = 1.0_dp c_E(7) = 1.0_dp ! A_E(2,1) = 1.0_dp/5.0_dp A_E(3,1) = 3.0_dp/40.0_dp A_E(3,2) = 9.0_dp/40.0_dp A_E(4,1) = 44.0_dp/45.0_dp A_E(4,2) = -56.0_dp/15.0_dp A_E(4,3) = 32.0_dp/9.0_dp A_E(5,1 ) = 19372.0_dp/6561.0_dp A_E(5,2) = -25360.0_dp/2187.0_dp A_E(5,3) = 64448.0_dp/6561.0_dp A_E(5,4) = -212.0_dp/729.0_dp A_E(6,1) = 9017.0_dp/3168.0_dp A_E(6,2)= -355.0_dp/33.0_dp A_E(6,3) = 46732.0_dp/5247.0_dp A_E(6,4) = 49.0_dp/176.0_dp A_E(6,5) = -5103.0_dp/18656.0_dp A_E(7,1) = 35.0_dp/384.0_dp A_E(7,3) = 500.0_dp/1113.0_dp A_E(7,4) = 125.0_dp/192.0_dp A_E(7,5) = -2187.0_dp/6784.0_dp A_E(7,6) = 11.0_dp/84.0_dp ! b_E(1) = 35.0_dp/384.0_dp b_E(2) = 0._dp b_E(3) = 500.0_dp/1113.0_dp b_E(4) = 125.0_dp/192.0_dp b_E(5) = -2187.0_dp/6784.0_dp b_E(6) = 11.0_dp/84.0_dp b_E(7) = 0._dp ! b_Es(1) = 5179.0_dp/57600.0_dp b_Es(2) = 0._dp b_Es(3) = 7571.0_dp/16695.0_dp b_Es(4) = 393.0_dp/640.0_dp b_Es(5) = -92097.0_dp/339200.0_dp b_Es(6) = 187.0_dp/2100.0_dp b_Es(7) = 1._dp/40._dp ! END SUBROUTINE DOPRI5_ADAPT END MODULE time_integration diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m index 8d494f0..e637d9d 100644 --- a/wk/analysis_gyacomo.m +++ b/wk/analysis_gyacomo.m @@ -1,233 +1,233 @@ %% UNCOMMENT FOR TUTORIAL gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory % resdir = '.'; %Name of the directory where the results are located % JOBNUMMIN = 00; JOBNUMMAX = 10; %% addpath(genpath([gyacomodir,'matlab'])) % ... add addpath(genpath([gyacomodir,'matlab/plot'])) % ... add addpath(genpath([gyacomodir,'matlab/compute'])) % ... add addpath(genpath([gyacomodir,'matlab/load'])) % ... add %% Load the results LOCALDIR = [gyacomodir,resdir,'/']; -MISCDIR = ['/misc/gyacomo_outputs/',resdir,'/']; %For long term storage +MISCDIR = ['/misc/',resdir,'/']; %For long term storage system(['mkdir -p ',MISCDIR]); system(['mkdir -p ',LOCALDIR]); CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD); system(CMD); % Load outputs from jobnummin up to jobnummax data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing data.localdir = LOCALDIR; data.FIGDIR = LOCALDIR; %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 1 %% Space time diagramm (fig 11 Ivanov 2020) % data.scale = 1;%/(data.Nx*data.Ny)^2; i_ = 1; disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))]) disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))]) options.TAVG_0 = data.TJOB_SE(i_)+600;%0.4*data.Ts3D(end); options.TAVG_1 = data.TJOB_SE(i_+1);%0.9*data.Ts3D(end); % Averaging times duration options.NCUT = 4; % Number of cuts for averaging and error estimation options.NMVA = 100; % Moving average for time traces % options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x) options.INTERP = 0; options.RESOLUTION = 256; fig = plot_radial_transport_and_spacetime(data,options); % save_figure(data,fig,'.png') end if 0 %% statistical transport averaging gavg =[]; gstd = []; for i_ = 3:2:numel(data.TJOB_SE) % i_ = 3; disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))]) disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))]) options.T = [data.TJOB_SE(i_) data.TJOB_SE(i_+1)]; options.NPLOTS = 0; [fig, gavg_, gstd_] = statistical_transport_averaging(data,options); gavg = [gavg gavg_]; gstd = [gstd gstd_]; end disp(gavg); disp(gstd); end if 0 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Options options.INTERP = 1; options.POLARPLOT = 0; -options.NAME = '\phi'; +% options.NAME = '\phi'; % options.NAME = '\omega_z'; -% options.NAME = 'N_i^{00}'; +options.NAME = 'N_e^{00}'; % options.NAME = 'v_x'; % options.NAME = 'n_i^{NZ}'; % options.NAME = '\Gamma_x'; % options.NAME = 'n_i'; -options.PLAN = 'kxky'; +options.PLAN = 'xy'; % options.NAME = 'f_i'; % options.PLAN = 'sx'; options.COMP = 'avg'; % options.TIME = data.Ts5D(end-30:end); % options.TIME = data.Ts3D; options.TIME = [000:1:8000]; data.EPS = 0.1; data.a = data.EPS * 2000; options.RESOLUTION = 256; create_film(data,options,'.gif') end if 0 %% 2D snapshots % Options options.INTERP = 1; options.POLARPLOT = 0; options.AXISEQUAL = 1; options.NAME = '\phi'; % options.NAME = '\psi'; % options.NAME = 'n_i'; % options.NAME = 'N_i^{00}'; % options.NAME = 'T_i'; % options.NAME = '\Gamma_x'; % options.NAME = 'k^2n_e'; options.PLAN = 'xy'; % options.NAME 'f_i'; % options.PLAN = 'sx'; options.COMP = 'avg'; options.TIME = [1000 3000 5000 8000 10000]; data.a = data.EPS * 2e3; fig = photomaton(data,options); % save_figure(data,fig) end if 0 %% 3D plot on the geometry options.INTERP = 0; options.NAME = '\phi'; options.PLANES = [1]; options.TIME = [30]; options.PLT_MTOPO = 1; options.PLT_FTUBE = 0; data.EPS = 0.4; data.rho_o_R = 3e-3; % Sound larmor radius over Machine size ratio fig = show_geometry(data,options); save_figure(data,fig,'.png') end if 0 %% Kinetic distribution function sqrt(xy) (GENE vsp) % options.SPAR = linspace(-3,3,32)+(6/127/2); % options.XPERP = linspace( 0,6,32); options.SPAR = gene_data.vp'; options.XPERP = gene_data.mu'; options.iz = 'avg'; options.T = [1500 5000]; options.PLT_FCT = 'contour'; options.ONED = 1; options.non_adiab = 0; options.SPECIE = 'e'; options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene fig = plot_fa(data,options); % save_figure(data,fig,'.png') end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 0; options.ST = 1; options.PLOT_TYPE = 'space-time'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [1000]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); % save_figure(data,fig,'.png'); end if 0 %% Time averaged spectrum options.TIME = [2000 3000]; options.NORM =1; % options.NAME = '\phi'; % options.NAME = 'N_i^{00}'; options.NAME ='\Gamma_x'; options.PLAN = 'kxky'; options.COMPZ = 'avg'; options.OK = 0; options.COMPXY = 'avg'; % avg/sum/max/zero/ 2D plot otherwise options.COMPT = 'avg'; options.PLOT = 'semilogy'; fig = spectrum_1D(data,options); % save_figure(data,fig,'.png') end if 0 %% 1D real plot options.TIME = [50 100 200]; options.NORM = 0; options.NAME = '\phi'; % options.NAME = 'n_i'; % options.NAME ='\Gamma_x'; % options.NAME ='s_y'; options.COMPX = 1; options.COMPY = 2; options.COMPZ = 1; options.COMPT = 1; options.MOVMT = 1; fig = real_plot_1D(data,options); % save_figure(data,fig,'.png') end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = [0.1 0.2 0.3 0.4]; options.TIME = [00:1200]; options.NMA = 1; options.NMODES = 5; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end if 0 %% ZF caracteristics (space time diagrams) TAVG_0 = 1200; TAVG_1 = 1500; % Averaging times duration % chose your field to plot in spacetime diag (uzf,szf,Gx) fig = ZF_spacetime(data,TAVG_0,TAVG_1); save_figure(data,fig,'.png') end if 0 %% Metric infos options.SHOW_FLUXSURF = 1; options.SHOW_METRICS = 0; fig = plot_metric(data,options); end if 0 %% linear growth rate for 3D fluxtube trange = [0 100]; nplots = 1; lg = compute_fluxtube_growth_rate(data,trange,nplots); end if 0 %% linear growth rate for 3D Zpinch trange = [5 20]; options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 1; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); % save_figure(data,fig,'.png') end diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m index 8a872de..5e9a28c 100644 --- a/wk/header_3D_results.m +++ b/wk/header_3D_results.m @@ -1,66 +1,71 @@ % Directory of the code "mypathtoHeLaZ/HeLaZ/" gyacomodir = '/home/ahoffman/gyacomo/'; % Directory of the simulation (from results) % if 1% Local results % resdir ='volcokas/64x32x16x5x3_kin_e_npol_1'; %% Dimits % resdir ='shearless_cyclone/128x64x16x5x3_Dim_90'; % resdir ='shearless_cyclone/128x64x16x9x5_Dim_scan/128x64x16x9x5_Dim_60'; % resdir ='shearless_cyclone/128x64x16x5x3_Dim_scan/128x64x16x5x3_Dim_70'; % resdir ='shearless_cyclone/64x32x16x5x3_Dim_scan/64x32x16x5x3_Dim_70'; %% AVS % resdir = 'volcokas/64x32x16x5x3_kin_e_npol_1'; % resdir = 'volcokas/64x32x16x5x3_kin_e_npol_1'; % resdir = 'shearless_cyclone/64x32x80x5x3_CBC_Npol_5_kine'; % resdir = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; % resdir = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine'; % resdir = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine'; %% shearless CBC % resdir ='shearless_cyclone/64x32x16x5x3_CBC_080'; % resdir ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100'; % resdir ='shearless_cyclone/64x32x16x5x3_CBC_120'; % resdir ='shearless_cyclone/64x32x16x9x5_CBC_080'; % resdir ='shearless_cyclone/64x32x16x9x5_CBC_100'; % resdir ='shearless_cyclone/64x32x16x9x5_CBC_120'; % resdir = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK'; %% CBC % resdir = 'CBC/64x32x16x5x3'; % resdir = 'CBC/64x128x16x5x3'; % resdir = 'CBC/128x64x16x5x3'; % resdir = 'CBC/96x96x16x3x2_Nexc_6'; % resdir = 'CBC/128x96x16x3x2'; % resdir = 'CBC/192x96x16x3x2'; -resdir = 'CBC/192x96x24x13x7'; +% resdir = 'CBC/192x96x24x13x7'; % resdir = 'CBC/kT_11_128x64x16x5x3'; % resdir = 'CBC/kT_9_256x128x16x3x2'; % resdir = 'CBC/kT_4.5_128x64x16x13x3'; % resdir = 'CBC/kT_4.5_192x96x24x13x7'; % resdir = 'CBC/kT_4.5_128x64x16x13x7'; % resdir = 'CBC/kT_4.5_128x96x24x15x5'; % resdir = 'CBC/kT_5.3_192x96x24x13x7'; % resdir = 'CBC/kT_13_large_box_128x64x16x5x3'; % resdir = 'CBC/kT_11_96x64x16x5x3_ky_0.02'; % resdir = 'CBC/kT_scan_128x64x16x5x3'; % resdir = 'CBC/kT_scan_192x96x16x3x2'; % resdir = 'CBC/kT_13_96x96x16x3x2_Nexc_6'; % resdir = 'dbg/nexc_dbg'; % resdir = 'CBC/NM_F4_kT_4.5_192x64x24x6x4'; % resdir = 'CBC_Ke_EM/192x96x24x5x3'; % resdir = 'CBC_Ke_EM/96x48x16x5x3'; % resdir = 'CBC_Ke_EM/minimal_res'; %% KBM % resdir = 'NL_KBM/192x64x24x5x3'; %% Linear CBC % resdir = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe'; % resdir = 'testcases/miller_example'; % resdir = 'Miller/128x256x3x2_CBC_dt_5e-3'; +%% CBC Miller +% resdir = 'GCM_CBC/daint/Miller_GX_comparison'; +%% RK3 +resdir = 'dbg/SSPx_RK3_test'; +resdir = 'dbg/SSPx_RK3_test/RK4'; resdir = ['results/',resdir]; JOBNUMMIN = 00; JOBNUMMAX = 10; run analysis_gyacomo diff --git a/wk/lin_ETPY.m b/wk/lin_ETPY.m index 101f81c..9b6c383 100644 --- a/wk/lin_ETPY.m +++ b/wk/lin_ETPY.m @@ -1,186 +1,187 @@ %% QUICK RUN SCRIPT for linear entropy mode (ETPY) in a Zpinch % This script create a directory in /results and run a simulation directly % from matlab framework. It is meant to run only small problems in linear % for benchmark and debugging purpose since it makes matlab "busy" % -SIMID = 'lin_ETPY'; % Name of the simulation -% SIMID = 'dbg'; % Name of the simulation -RUN = 0; % To run or just to load +% SIMID = 'lin_ETPY'; % Name of the simulation +SIMID = 'dbg'; % Name of the simulation +RUN = 1; % To run or just to load addpath(genpath('../matlab')) % ... add default_plots_options PROGDIR = '/home/ahoffman/gyacomo/'; % EXECNAME = 'gyacomo_dbg'; EXECNAME = 'gyacomo'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS -NU = 0.01; % Collision frequency +NU = 0.05; % Collision frequency TAU = 1.0; % e/i temperature ratio -K_Ne = 2.5; % ele Density ''' -K_Te = K_Ne/4; % ele Temperature ''' -K_Ni = K_Ne; % ion Density gradient drive -K_Ti = K_Ni/4; % ion Temperature ''' +K_Ne = 2.0; % ele Density ''' +K_Te = 0.4; % ele Temperature ''' +K_Ni = 2.0; % ion Density gradient drive +K_Ti = 0.4; % ion Temperature ''' SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons BETA = 0.0; % electron plasma beta %% GRID PARAMETERS P = 4; J = P/2; PMAXE = P; % Hermite basis size of electrons JMAXE = J; % Laguerre " PMAXI = P; % " ions JMAXI = J; % " NX = 2; % real space x-gridpoints -NY = 100; % '' y-gridpoints +NY = 50; % '' y-gridpoints LX = 2*pi/0.8; % Size of the squared frequency domain -LY = 200;%2*pi/0.05; % Size of the squared frequency domain +LY = 2*pi/0.05; % Size of the squared frequency domain NZ = 1; % number of perpendicular planes (parallel grid) NPOL = 1; SG = 0; % Staggered z grids option %% GEOMETRY -GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps +GEOMETRY= 'Z-pinch'; Q0 = 1.4; % safety factor SHEAR = 0.8; % magnetic shear -KAPPA = 0.0; % elongation +KAPPA = 1.0; % elongation DELTA = 0.0; % triangularity ZETA = 0.0; % squareness NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) EPS = 0.18; % inverse aspect ratio %% TIME PARMETERS -TMAX = 500; % Maximal time unit -DT = 5e-2; % Time step +TMAX = 100; % Maximal time unit +DT = 1e-2; % Time step SPS0D = 1; % Sampling per time unit for 2D arrays SPS2D = 0; % Sampling per time unit for 2D arrays SPS3D = 1; % Sampling per time unit for 2D arrays SPS5D = 1/5; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints JOB2LOAD= -1; %% OPTIONS LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) % Collision operator % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) -CO = 'LD'; -GKCO = 1; % gyrokinetic operator +CO = 'DG'; +GKCO = 0; % gyrokinetic operator ABCO = 1; % interspecies collisions INIT_ZF = 0; ZF_AMP = 0.0; CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS) KERN = 0; % Kernel model (0 : GK) -INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom +INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom %% OUTPUTS W_DOUBLE = 1; W_GAMMA = 1; W_HF = 1; W_PHI = 1; W_NA00 = 1; W_DENS = 1; W_TEMP = 1; W_NAPJ = 1; W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % unused HD_CO = 0.0; % Hyper diffusivity cutoff ratio MU = 0.0; % Hyperdiffusivity coefficient INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; MU_X = MU; % MU_Y = MU; % N_HD = 4; -MU_Z = 2.0; % +MU_Z = 0.0; % MU_P = 0.0; % MU_J = 0.0; % LAMBDAD = 0.0; NOISE0 = 1.0e-5; % Init noise amplitude BCKGD0 = 0.0; % Init background GRADB = 1.0; CURVB = 1.0; COLL_KCUT = 1000; %%------------------------------------------------------------------------- %% RUN setup % system(['rm fort*.90']); % Run linear simulation if RUN system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',PROGDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk']) end %% Load results %% filename = [SIMID,'/',PARAMS,'/']; LOCALDIR = [PROGDIR,'results/',filename,'/']; % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 00; data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing %% Short analysis if 1 %% linear growth rate (adapted for 2D zpinch and fluxtube) options.TRANGE = [0.5 1]*data.Ts3D(end); options.NPLOTS = 2; % 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 lg = compute_fluxtube_growth_rate(data,options); [gmax, kmax] = max(lg.g_ky(:,end)); [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky); msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg); msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg); end if 0 %% Ballooning plot options.time_2_plot = [120]; options.kymodes = [0.9]; options.normalized = 1; % options.field = 'phi'; fig = plot_ballooning(data,options); end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 1; options.ST = 1; options.PLOT_TYPE = 'space-time'; % options.PLOT_TYPE = 'Tavg-1D'; % options.PLOT_TYPE = 'Tavg-2D'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [0 50]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); save_figure(data,fig) end if 0 %% linear growth rate for 3D Zpinch (kz fourier transform) trange = [0.5 1]*data.Ts3D(end); options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 0; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig) end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = 10; options.TIME = [0:1000]; options.NMA = 1; options.NMODES = 10; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end if 0 %% RH TEST ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end); [~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D)); plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3)); figure plot(data.Ts3D(it0:it1), plt(data.PHI)); xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$') title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky))) end diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m index 6f18bee..ebca810 100644 --- a/wk/lin_ITG.m +++ b/wk/lin_ITG.m @@ -1,190 +1,191 @@ %% QUICK RUN SCRIPT % This script create a directory in /results and run a simulation directly % from matlab framework. It is meant to run only small problems in linear % for benchmark and debugging purpose since it makes matlab "busy" % -SIMID = 'lin_ITG'; % Name of the simulation -RUN = 0; % To run or just to load +% SIMID = 'lin_ITG'; % Name of the simulation +SIMID = 'dbg'; % Name of the simulation +RUN = 1; % To run or just to load addpath(genpath('../matlab')) % ... add default_plots_options gyacomodir = '/home/ahoffman/gyacomo/'; % EXECNAME = 'gyacomo_1.0'; % EXECNAME = 'gyacomo_dbg'; EXECNAME = 'gyacomo'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS NU = 0.05; % Collision frequency TAU = 1.0; % e/i temperature ratio -K_Ne = 2.22; % ele Density ''' +K_Ne = 0*2.22; % ele Density ''' K_Te = 6.96; % ele Temperature ''' -K_Ni = 2.22; % ion Density gradient drive +K_Ni = 0*2.22; % ion Density gradient drive K_Ti = 6.96; % ion Temperature ''' SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) % SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380) KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons BETA = 0.0; % electron plasma beta %% GRID PARAMETERS P = 4; J = P/2; PMAXE = P; % Hermite basis size of electrons JMAXE = J; % Laguerre " PMAXI = P; % " ions JMAXI = J; % " NX = 20; % real space x-gridpoints NY = 2; % '' y-gridpoints LX = 2*pi/0.8; % Size of the squared frequency domain LY = 2*pi/0.3; % Size of the squared frequency domain NZ = 32; % number of perpendicular planes (parallel grid) NPOL = 1; SG = 0; % Staggered z grids option %% GEOMETRY % GEOMETRY= 's-alpha'; GEOMETRY= 'miller'; Q0 = 1.4; % safety factor SHEAR = 0.8; % magnetic shear -KAPPA = 0.0; % elongation +KAPPA = 1.0; % elongation DELTA = 0.0; % triangularity ZETA = 0.0; % squareness NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear)) EPS = 0.18; % inverse aspect ratio %% TIME PARMETERS -TMAX = 1; % Maximal time unit +TMAX = 25; % Maximal time unit DT = 3e-3; % Time step SPS0D = 1; % Sampling per time unit for 2D arrays SPS2D = 0; % Sampling per time unit for 2D arrays SPS3D = 5; % Sampling per time unit for 2D arrays SPS5D = 1/5; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints JOB2LOAD= -1; %% OPTIONS LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1) % Collision operator % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau) CO = 'DG'; GKCO = 0; % gyrokinetic operator ABCO = 1; % interspecies collisions INIT_ZF = 0; ZF_AMP = 0.0; CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS) KERN = 0; % Kernel model (0 : GK) INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom %% OUTPUTS W_DOUBLE = 1; W_GAMMA = 1; W_HF = 1; W_PHI = 1; W_NA00 = 1; W_DENS = 1; W_TEMP = 1; W_NAPJ = 1; W_SAPJ = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % unused HD_CO = 0.0; % Hyper diffusivity cutoff ratio MU = 0.0; % Hyperdiffusivity coefficient INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; MU_X = MU; % MU_Y = MU; % N_HD = 4; MU_Z = 2.0; % MU_P = 0.0; % MU_J = 0.0; % LAMBDAD = 0.0; NOISE0 = 0.0e-5; % Init noise amplitude BCKGD0 = 1.0; % Init background GRADB = 1.0; CURVB = 1.0; %%------------------------------------------------------------------------- %% RUN setup % system(['rm fort*.90']); % Run linear simulation if RUN % system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk']) system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk']) end %% Load results %% filename = [SIMID,'/',PARAMS,'/']; LOCALDIR = [gyacomodir,'results/',filename,'/']; % Load outputs from jobnummin up to jobnummax JOBNUMMIN = 00; JOBNUMMAX = 00; data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing %% Short analysis if 1 %% linear growth rate (adapted for 2D zpinch and fluxtube) options.TRANGE = [0.5 1]*data.Ts3D(end); options.NPLOTS = 2; % 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 lg = compute_fluxtube_growth_rate(data,options); [gmax, kmax] = max(lg.g_ky(:,end)); [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky); msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg); msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg); end if 1 %% Ballooning plot options.time_2_plot = [120]; options.kymodes = [0.3]; options.normalized = 1; % options.field = 'phi'; fig = plot_ballooning(data,options); end if 0 %% Hermite-Laguerre spectrum % options.TIME = 'avg'; options.P2J = 1; options.ST = 1; options.PLOT_TYPE = 'space-time'; % options.PLOT_TYPE = 'Tavg-1D'; % options.PLOT_TYPE = 'Tavg-2D'; options.NORMALIZED = 0; options.JOBNUM = 0; options.TIME = [0 50]; options.specie = 'i'; options.compz = 'avg'; fig = show_moments_spectrum(data,options); % fig = show_napjz(data,options); save_figure(data,fig) end if 0 %% linear growth rate for 3D Zpinch (kz fourier transform) trange = [0.5 1]*data.Ts3D(end); options.keq0 = 1; % chose to plot planes at k=0 or max options.kxky = 1; options.kzkx = 0; options.kzky = 0; [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options); save_figure(data,fig) end if 0 %% Mode evolution options.NORMALIZED = 0; options.K2PLOT = 1; options.TIME = [0:1000]; options.NMA = 1; options.NMODES = 1; options.iz = 'avg'; fig = mode_growth_meter(data,options); save_figure(data,fig,'.png') end if 1 %% RH TEST ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end); [~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D)); plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3)); figure plot(data.Ts3D(it0:it1), plt(data.PHI)); xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$') title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky))) end