Add triad energy transfer diagnostics

- The code does not implement OPENACC or OMPGPU. - Set triad_print_flag=1 to save the bin.cgyro.triad file. - Saved data includes : entropy, field energy, triad energy transfer (T), triad energy transfer with only non-zonal pairs (T'), and dissipation (radial, theta, collisions). - Visualization options : -plot triad, triad_v2
gafusion · Nov 4, 2024 · dfc8b32 · dfc8b32
1 parent 4114b39
commit dfc8b32
Show file tree

Hide file tree

Showing 10 changed files with 785 additions and 5 deletions.
diff --git a/cgyro/src/cgyro_flux.f90 b/cgyro/src/cgyro_flux.f90
@@ -44,6 +44,7 @@ subroutine cgyro_flux
   complex :: cprod
   real, parameter :: x_fraction=0.2
   real :: u
+  real :: kx
 
 !$omp parallel do private(iv_loc,iv,is,ix,ie,dv,vpar,ic,ir,it,erot,cprod,cn) &
 !$omp&            private(prod1,prod2,prod3,l,icl,dvr,u,flux_norm) &
@@ -188,6 +189,34 @@ subroutine cgyro_flux
         cflux_loc(:,:,:,itor) = cflux_loc(:,:,:,itor)+2*sin(u)*real(gflux_loc(l,:,:,:,itor))/u
      enddo
 
+
+  !-------------------------------------------------------------
+  ! 2-1. Compute Triad energy transfer
+  !-------------------------------------------------------------
+
+     kx = 2*pi*rho/length
+     do is=1,n_species
+        ! Triad energy transfer   :   T_k
+        triad_loc_old(is,:,itor,1) = triad_loc(is,:,itor,1) *temp(is)/dlntdr(is_ele)
+        ! From Nonzonal Triad energy transfer   :   T_k [NZ(k',k")->k]
+        triad_loc_old(is,:,itor,2) = triad_loc(is,:,itor,2) *temp(is)/dlntdr(is_ele)
+        ! dEntropy / dt   :   dS_k/dt 
+        triad_loc_old(is,:,itor,3) = (triad_loc(is,:,itor,3)-triad_loc_old(is,:,itor,3))/delta_t
+        triad_loc_old(is,:,itor,3) = triad_loc_old(is,:,itor,3) *temp(is)/dlntdr(is_ele)*0.5
+        ! dWk_perp / dt    , (ky = 0)  
+        triad_loc_old(is,:,itor,4) = (triad_loc(is,:,itor,4)-triad_loc_old(is,:,itor,4))/delta_t
+        triad_loc_old(is,:,itor,4) = triad_loc_old(is,:,itor,4) *(kx*lambda_star)**2/dlntdr(is_ele)*0.5
+        ! Entropy   :   S_k
+        triad_loc_old(is,:,itor,5) = triad_loc(is,:,itor,3) *temp(is)/dlntdr(is_ele)*0.5
+        ! Diss. (radial)
+        triad_loc_old(is,:,itor,6) = triad_loc(is,:,itor,5) *temp(is)/dlntdr(is_ele)
+        ! Diss. (theta )
+        triad_loc_old(is,:,itor,7) = triad_loc(is,:,itor,6) *temp(is)/dlntdr(is_ele)
+        ! Diss. (Coll. )
+        triad_loc_old(is,:,itor,8) = triad_loc(is,:,itor,7) *temp(is)/dlntdr(is_ele)
+     enddo
+
+
      !-----------------------------------------------------
      ! 3. Renormalize fluxes to GB or quasilinear forms
      !~----------------------------------------------------
@@ -216,7 +245,7 @@ subroutine cgyro_flux
         ! GyroBohm normalizations
         gflux_loc(:,:,:,:,itor) = gflux_loc(:,:,:,:,itor)/rho**2
         cflux_loc(:,:,:,itor) = cflux_loc(:,:,:,itor)/rho**2
-
+        triad_loc_old(:,:,itor,:)  = triad_loc_old(:,:,itor,:)/rho**2
      endif
   enddo
 
@@ -240,6 +269,15 @@ subroutine cgyro_flux
        NEW_COMM_1, &
        i_err)
 
+  ! Reduced complex triad(ns,kx), below, is still distributed over n 
+  call MPI_ALLREDUCE(triad_loc_old(:,:,:,:), &
+       triad, &
+       size(triad), &
+       MPI_DOUBLE_COMPLEX, &
+       MPI_SUM, &
+       NEW_COMM_1, &
+       i_err)
+
   ! Reduced real cflux(ky), below, is still distributed over n 
 
   call MPI_ALLREDUCE(cflux_loc, &

diff --git a/cgyro/src/cgyro_globals.F90 b/cgyro/src/cgyro_globals.F90
@@ -80,6 +80,7 @@ module cgyro_globals
   integer :: moment_print_flag
   integer :: gflux_print_flag
   integer :: field_print_flag
+  integer :: triad_print_flag
   real :: amp0
   real :: amp
   real :: gamma_e
@@ -192,7 +193,7 @@ module cgyro_globals
   integer :: nsplit,nsplitA,nsplitB
   integer :: ns1,ns2
   integer, dimension(:), allocatable :: recv_status
-  integer :: fA_req, fB_req, g_req
+  integer :: fA_req, fB_req, eA_req, eB_req, g_req
   logical :: fA_req_valid, fB_req_valid, g_req_valid
   ! Thetas present in the process after NL AllToAll
   integer :: n_jtheta
@@ -243,6 +244,7 @@ module cgyro_globals
   character(len=12) :: binfile_hb      = 'bin.cgyro.hb'
   character(len=17) :: binfile_ky_flux = 'bin.cgyro.ky_flux'
   character(len=18) :: binfile_ky_cflux = 'bin.cgyro.ky_cflux'
+  character(len=15) :: binfile_triad   = 'bin.cgyro.triad'
   character(len=15), dimension(3) :: binfile_fieldb = &
        (/'bin.cgyro.phib ','bin.cgyro.aparb','bin.cgyro.bparb'/)
   character(len=16), dimension(3) :: binfile_kxky = &
@@ -344,11 +346,14 @@ module cgyro_globals
   complex, dimension(:,:,:), allocatable :: h0_x
   complex, dimension(:,:,:), allocatable :: h0_old
   complex, dimension(:,:,:,:), allocatable :: fA_nl,fB_nl
+  complex, dimension(:,:,:,:), allocatable :: eA_nl,eB_nl
   complex, dimension(:,:,:,:), allocatable :: g_nl
   complex, dimension(:,:,:), allocatable :: fpackA,fpackB
+  complex, dimension(:,:,:), allocatable :: epackA,epackB
   complex, dimension(:,:,:,:), allocatable :: gpack
   complex, dimension(:,:,:), allocatable :: omega_cap_h
   complex, dimension(:,:,:), allocatable :: omega_h
+  complex, dimension(:,:,:), allocatable :: diss_r
   complex, dimension(:,:,:,:), allocatable :: omega_s,omega_ss
   complex, dimension(:,:,:), allocatable :: omega_sbeta
   complex, dimension(:,:,:), allocatable :: cap_h_c
@@ -382,6 +387,7 @@ module cgyro_globals
   real, dimension(:,:,:,:), allocatable :: cflux
   complex, dimension(:,:,:,:,:), allocatable :: gflux_loc
   complex, dimension(:,:,:,:,:), allocatable :: gflux
+  complex, dimension(:,:,:,:), allocatable :: triad,triad_loc,triad_loc_old
   real, dimension(:,:), allocatable :: cflux_tave, gflux_tave
   real :: tave_min, tave_max
   integer :: tave_step

diff --git a/cgyro/src/cgyro_init_arrays.F90 b/cgyro/src/cgyro_init_arrays.F90
@@ -448,7 +448,14 @@ subroutine cgyro_init_arrays
                 + abs(omega_cdrift_r(it,is)*xi(ix))*vel(ie) &
                 + abs(omega_rot_drift_r(it,is)) &
                 + abs(omega_rot_edrift_r(it)))          
-
+
+            ! (d/dr) upwind diss.
+            diss_r(ic,iv_loc,itor) = - (n_radial/length)*spectraldiss(u,nup_radial)*up_radial &
+             * (abs(omega_rdrift(it,is))*energy(ie)*(1.0+xi(ix)**2) &
+             + abs(omega_cdrift_r(it,is)*xi(ix))*vel(ie) &
+             + abs(omega_rot_drift_r(it,is)) &
+             + abs(omega_rot_edrift_r(it))) 
+
            ! omega_star 
            carg = &
                 -i_c*k_theta_base*itor*rho*(dlnndr(is)+dlntdr(is)*(energy(ie)-1.5)) &

diff --git a/cgyro/src/cgyro_init_manager.F90 b/cgyro/src/cgyro_init_manager.F90
@@ -178,6 +178,11 @@ subroutine cgyro_init_manager
      allocate(cflux_loc(n_species,4,n_field,nt1:nt2))
      allocate(    gflux(0:n_global,n_species,4,n_field,nt1:nt2))
      allocate(gflux_loc(0:n_global,n_species,4,n_field,nt1:nt2))
+
+     allocate(    triad(n_species,n_radial,nt1:nt2,8))
+     allocate(triad_loc(n_species,n_radial,nt1:nt2,7))
+     allocate(triad_loc_old(n_species,n_radial,nt1:nt2,8))
+
      allocate(cflux_tave(n_species,4))
      allocate(gflux_tave(n_species,4))
 
@@ -245,6 +250,7 @@ subroutine cgyro_init_manager
      allocate(cap_h_v(nc_loc,nt1:nt2,nv))
      allocate(omega_cap_h(nc,nv_loc,nt1:nt2))
      allocate(omega_h(nc,nv_loc,nt1:nt2))
+     allocate(diss_r(nc,nv_loc,nt1:nt2))
      allocate(omega_s(n_field,nc,nv_loc,nt1:nt2))
      allocate(omega_ss(n_field,nc,nv_loc,nt1:nt2))
      allocate(omega_sbeta(nc,nv_loc,nt1:nt2))
@@ -285,8 +291,10 @@ subroutine cgyro_init_manager
      ! Nonlinear arrays
      if (nonlinear_flag == 1) then
         allocate(fA_nl(n_radial,nt_loc,nsplitA,n_toroidal_procs))
+        allocate(eA_nl(n_radial,nt_loc,nsplitA,n_toroidal_procs))
         allocate(g_nl(n_field,n_radial,n_jtheta,n_toroidal))
         allocate(fpackA(n_radial,nt_loc,nsplitA*n_toroidal_procs))
+        allocate(epackA(n_radial,nt_loc,nsplitA*n_toroidal_procs))
         allocate(gpack(n_field,n_radial,n_jtheta,n_toroidal))
         allocate(jvec_c_nl(n_field,n_radial,n_jtheta,nv_loc,n_toroidal))
 #if defined(OMPGPU)
@@ -297,6 +305,8 @@ subroutine cgyro_init_manager
         if (nsplitB > 0) then ! nsplitB can be zero at large MPI
           allocate(fB_nl(n_radial,nt_loc,nsplitB,n_toroidal_procs))
           allocate(fpackB(n_radial,nt_loc,nsplitB*n_toroidal_procs))
+          allocate(eB_nl(n_radial,nt_loc,nsplitB,n_toroidal_procs))
+          allocate(epackB(n_radial,nt_loc,nsplitB*n_toroidal_procs))
 #if defined(OMPGPU)
 !$omp target enter data map(alloc:fpackB,fB_nl)
 #elif defined(_OPENACC)