collisional_heating.f90 Source File


Source Code

!> FIXME : Add documentation
module collisional_heating

  implicit none
  public :: init_collisional, finish_collisional
  public :: write_collisional, calculate_collisional
  complex, dimension (:,:,:), allocatable :: g0   !< g0 will be used in calculate_collisional,
  complex, dimension (:,:,:,:), allocatable :: tot              !< they must be dummy arguments...
  complex, dimension(:,:), allocatable :: coll_heating, coll_heating_2 !< we will use the hk variable in diagnostics_heating
  logical :: initialized = .false.

  !> FIXME : Add documentation  
  subroutine init_collisional(gnostics) 
    use kt_grids, only: naky, ntheta0, nx, ny
    use theta_grid, only: ntgrid
    use species, only: nspec
    use diagnostics_config, only: diagnostics_type
    use le_grids, only: init_le_grids, nlambda, negrid
    use gs2_transforms, only: init_transforms
    use gs2_layouts, only: g_lo, init_dist_fn_layouts
    use hyper, only: init_hyper
    use mp, only: nproc, iproc
    implicit none
    type(diagnostics_type), intent(inout) :: gnostics
    logical :: accel_x

    if(.not. gnostics%write_collisional .or. initialized) return
    initialized = .true.

!!!calls every initiation routine we need: le_grids for integration, init_gs2 transforms to initiate fourier transforms routines...
    ! benefit: calls init_theta_grid, init_kt_grids, init_gs2_layouts and init_species for later use
    call init_le_grids
    ! in order to use fft transform routines
    call init_transforms(ntgrid, naky, ntheta0, nlambda, negrid, nspec, nx, ny, accel_x)
    call init_dist_fn_layouts(ntgrid, naky, ntheta0, nlambda, negrid, nspec, nproc, iproc)

!!!!! now need to initiate and assign hypervisc_filter
    call init_hyper
!!!!! allocate and initialize arrays
    allocate (coll_heating(ntheta0,naky))
    allocate (coll_heating_2(ntheta0,naky))
  end subroutine init_collisional

  !> FIXME : Add documentation  
  subroutine finish_collisional
    if(.not. initialized) return
    initialized = .false.
  end subroutine finish_collisional

  !> FIXME : Add documentation  
  subroutine calculate_collisional()
    use dist_fn_arrays, only: g, gnew
    use diagnostics_config, only: diagnostics_type
    use dist_fn_arrays, only: c_rate, g_adjust, from_g_gs2
    use theta_grid, only: ntgrid, delthet, jacob
    use kt_grids, only: naky, ntheta0, aky
    use species, only: spec, nspec
    use le_grids, only: integrate_moment
    use gs2_layouts, only: g_lo, it_idx, ik_idx, is_idx
    use gs2_time, only: code_dt
    use hyper, only: hyper_diff, hypervisc_filter
    use fields_arrays, only: phi, bpar
    use collisions, only: heating

    implicit none

    real, dimension(-ntgrid:ntgrid) :: wgt !< weights for theta integration
    complex, dimension(ntheta0,naky) :: hyper_part, hyper_part_2, colls_part   !< two parts of c_rate array: collisions and hypervisocity
    !complex, dimension(-ntgrid:ntgrid,ntheta0,naky) :: tot
    integer :: is, ik, it, ig, isgn, iglo
    complex :: dgdt_hypervisc, havg
    real :: fac, fac2
    real, dimension(1:nspec) :: weights 

    hyper_part = 0.
    hyper_part_2 = 0.
    colls_part = 0.

    wgt = delthet*jacob
    wgt = wgt/sum(wgt)
    call g_adjust(g0,phi,bpar, direction = from_g_gs2)   !< transforms g into h (minus boltzmann response), anyway fbpar=0

    call hyper_diff(g0,phi) !< assigns values to hypervisc_filter
    if (heating .and. allocated(c_rate)) then
       do is = 1, nspec
          do ik = 1, naky
             fac = 0.5
             if (aky(ik) < epsilon(0.)) fac = 1.0
             do it = 1, ntheta0
                do ig= -ntgrid,ntgrid
                   colls_part(it,ik) = colls_part(it,ik) + real(c_rate(ig,it,ik,is,3))*fac*wgt(ig)*spec(is)%temp*spec(is)%dens
                   hyper_part_2(it,ik) = hyper_part_2(it,ik) + real(c_rate(ig,it,ik,is,2))*fac*wgt(ig)*spec(is)%temp*spec(is)%dens
                end do
             end do
          end do
       end do
    end if

    do iglo=g_lo%llim_proc, g_lo%ulim_proc
      is = is_idx(g_lo, iglo)
      it = it_idx(g_lo, iglo)
      ik = ik_idx(g_lo, iglo)
      do isgn=1,2
        do ig=-ntgrid, ntgrid-1
          havg = 0.25*(g(ig,isgn,iglo)+ g(ig+1,isgn,iglo)+ gnew(ig,isgn,iglo)+ gnew(ig+1,isgn,iglo)) 

          dgdt_hypervisc = 0.5*((1.0-1./hypervisc_filter(ig,it,ik))*gnew(ig,isgn,iglo) &
                     + (1.0-1./hypervisc_filter(ig+1,it,ik))*gnew(ig+1,isgn,iglo))/code_dt

          g0(ig,isgn,iglo) = spec(is)%dens*spec(is)%temp*conjg(havg)*dgdt_hypervisc
        end do
      end do
    end do

    !call integrate_species_original(g0, weights, tot)    !!!!! integrate species_original is dubious...
    call integrate_moment(g0, tot)
    !write(*,*) tot
    do ik = 1, naky
       fac2 = 0.5
       if (aky(ik) < epsilon(0.0)) fac2 = 1.0
       do it = 1, ntheta0
          do is = 1, nspec
             do ig = -ntgrid, ntgrid-1
                hyper_part(it,ik) = hyper_part(it,ik) + real(tot(ig,it,ik,is))*wgt(ig)*fac2
             end do
          end do
       end do
    end do
    coll_heating= colls_part + hyper_part
    coll_heating_2 = colls_part + hyper_part_2
  end subroutine calculate_collisional

  !> FIXME : Add documentation  
  subroutine write_collisional(gnostics)
    use gs2_io, only: starts, netcdf_write_complex, omega_dims
    use fields_parallelization, only: field_k_local
    use diagnostics_config, only: diagnostics_type
    use mp, only: proc0
    implicit none
    type(diagnostics_type), intent(in) :: gnostics
    if (.not. proc0) return

    call netcdf_write_complex(gnostics%file_id, "collisional_heating", coll_heating, &
         dim_names=omega_dims, start=starts(4, gnostics%nout),  &
         long_name="Collisional and hyper viscous rate of loss of free energy for each mode and each t")
    call netcdf_write_complex(gnostics%file_id, "collisional_heating_2", coll_heating_2, &
         dim_names=omega_dims, start=starts(4, gnostics%nout),  &
         long_name="Collisional and hyper viscous rate of loss of free energy for each mode and each t, other method")
  end subroutine write_collisional

end module collisional_heating