!> FIXME : Add documentation module collisional_heating implicit none private 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. contains !> 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(tot(-ntgrid:ntgrid,ntheta0,naky,nspec)) allocate(g0(-ntgrid:ntgrid,2,g_lo%llim_proc:g_lo%ulim_alloc)) allocate (coll_heating(ntheta0,naky)) coll_heating=0.0 allocate (coll_heating_2(ntheta0,naky)) coll_heating_2=0.0 end subroutine init_collisional !> FIXME : Add documentation subroutine finish_collisional if(.not. initialized) return deallocate(coll_heating) deallocate(coll_heating_2) deallocate(tot) deallocate(g0) 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) g0=g weights=1. 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