Calculate the collision error and write to text file <runname>.cres
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| complex, | intent(in), | dimension (-ntgrid:,:,:) | :: | phi | ||
| complex, | intent(in), | dimension (-ntgrid:,:,:) | :: | bpar | ||
| logical, | intent(in) | :: | last |
If true, close the file |
subroutine collision_error (phi, bpar, last)
use mp, only: proc0, send, receive, barrier
use le_grids, only: ng2, jend, nlambda, lambda_map
use theta_grid, only: ntgrid
use dist_fn_arrays, only: gnew, g_adjust, g_work, to_g_gs2, from_g_gs2
use gs2_layouts, only: lz_lo, ig_idx, idx_local, proc_id
use gs2_layouts, only: ik_idx, ie_idx, is_idx, it_idx, il_idx
use collisions, only: dtot, fdf, fdb
use redistribute, only: gather, scatter
use file_utils, only: open_output_file, close_output_file
use gs2_time, only: user_time
use array_utils, only: zero_array
implicit none
complex, dimension (-ntgrid:,:,:), intent (in) :: phi, bpar
!> If true, close the file
logical, intent (in) :: last
integer :: je, te, ig, il, ip, ilz, ie, is, ik, it
integer :: igmax, ikmax, itmax, iemax, ilmax, ismax
integer, save :: unit
complex, dimension (:), allocatable :: ltmp, ftmp
complex, dimension (:,:), allocatable :: lcoll, fdcoll, glze
! logical :: first = .true.
real :: etmp, emax, etot, eavg, edenom, ltmax
real :: time
allocate (ltmp(2*nlambda), ftmp(2*nlambda))
allocate (lcoll(2*nlambda,lz_lo%llim_proc:lz_lo%ulim_alloc))
allocate (fdcoll(2*nlambda,lz_lo%llim_proc:lz_lo%ulim_alloc))
allocate (glze(max(2*nlambda,2*ng2+1),lz_lo%llim_proc:lz_lo%ulim_alloc))
call zero_array(lcoll); call zero_array(fdcoll); call zero_array(glze)
ltmp = 0.0; ftmp = 0.0
etmp = 0.0; emax = 0.0; etot = 0.0; eavg = 0.0; edenom = 1.0; ltmax = 1.0
! if (first .and. proc0) then
if (.not.cerrinit .and. proc0) then
call open_output_file (unit,".cres")
! first = .false.
cerrinit = .true.
end if
! convert gnew from g to h
call g_adjust(gnew,phi,bpar, direction = from_g_gs2)
g_work = gnew
! convert gnew from h back to g
call g_adjust(gnew,phi,bpar, direction = to_g_gs2)
! map from g_work(ig,isgn,iglo) to glze(il,ilz)
call gather (lambda_map, g_work, glze)
! loop over ig, isign, ik, it, ie, is
do ilz = lz_lo%llim_proc, lz_lo%ulim_proc
ig = ig_idx(lz_lo,ilz)
if (jend(ig) == 0) then ! no trapped particles
! je = number of + vpa grid pts
! te = number of + and - vpa grid pts
je = ng2
te = 2*ng2
! find d/d(xi) ((1+xi**2)( d g(xi)/ d(xi) )) at each xi
! using lagrange (lcoll) and finite difference (fdcoll)
il = 1
do ip = il, il+2
lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,il,ip)*glze(ip,ilz)
end do
il = 2
do ip = il-1, il+1
lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,il,ip-il+2)*glze(ip,ilz)
end do
do il=3,ng2
do ip=il-2,il+2
lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,il,ip-il+3)*glze(ip,ilz)
end do
end do
do il=ng2+1, 2*ng2-2
do ip = il-2,il+2
lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,2*ng2-il+1,il-ip+3)*glze(ip,ilz)
end do
end do
il = 2*ng2-1
do ip = il-1, il+1
lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,2,il-ip+2)*glze(ip,ilz)
end do
il = 2*ng2
do ip = il-2, il
lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,1,il-ip+1)*glze(ip,ilz)
end do
! deal with xi from 1-eps -> eps
do il=2,ng2
fdcoll(il,ilz) = fdf(ig,il)*(glze(il+1,ilz) - glze(il,ilz)) - fdb(ig,il)*(glze(il,ilz) - glze(il-1,ilz))
end do
! deal with xi from -eps -> -1+eps
do il=ng2+1, 2*ng2-1
fdcoll(il,ilz) = fdb(ig,2*ng2-il+1)*(glze(il+1,ilz) - glze(il,ilz)) - fdf(ig,2*ng2-il+1)*(glze(il,ilz) - glze(il-1,ilz))
end do
fdcoll(1,ilz) = fdf(ig,1)*(glze(2,ilz) - glze(1,ilz))
fdcoll(2*ng2,ilz) = -fdf(ig,1)*(glze(2*ng2,ilz) - glze(2*ng2-1,ilz))
else ! trapped particle runs
je = jend(ig)
te = 2*je - 1
il = 1
do ip = il, il+2
lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,il,ip)*glze(ip,ilz)
end do
il = 2
do ip = il-1, il+1
lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,il,ip-il+2)*glze(ip,ilz)
end do
do il=3,je
do ip=il-2,il+2
lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,il,ip-il+3)*glze(ip,ilz)
end do
end do
do il=je+1, te-2
do ip = il-2,il+2
lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,te-il+1,il-ip+3)*glze(ip,ilz)
end do
end do
il = te-1
do ip = il-1, il+1
lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,2,il-ip+2)*glze(ip,ilz)
end do
il = te
do ip = il-2, il
lcoll(il,ilz) = lcoll(il,ilz) + dtot(ig,1,il-ip+1)*glze(ip,ilz)
end do
! is il=je handled correctly here?
do il=2,je
fdcoll(il,ilz) = fdf(ig,il)*(glze(il+1,ilz) - glze(il,ilz)) - fdb(ig,il)*(glze(il,ilz) - glze(il-1,ilz))
end do
do il=je+1,te-1
fdcoll(il,ilz) = fdb(ig,te-il+1)*(glze(il+1,ilz) - glze(il,ilz)) - fdf(ig,te-il+1)*(glze(il,ilz) - glze(il-1,ilz))
end do
fdcoll(1,ilz) = fdf(ig,1)*(glze(2,ilz) - glze(1,ilz))
fdcoll(te,ilz) = -fdf(ig,1)*(glze(te,ilz) - glze(te-1,ilz))
end if
end do
time = user_time
do ilz=lz_lo%llim_world, lz_lo%ulim_world
ig = ig_idx(lz_lo, ilz)
ik = ik_idx(lz_lo, ilz)
it = it_idx(lz_lo, ilz)
ie = ie_idx(lz_lo, ilz)
is = is_idx(lz_lo, ilz)
je = jend(ig)
if (je == 0) then
te = 2*ng2
else
te = 2*je-1
end if
if (idx_local (lz_lo, ilz)) then
if (proc0) then
ltmp = lcoll(:,ilz)
ftmp = fdcoll(:,ilz)
else
call send (lcoll(:,ilz), 0)
call send (fdcoll(:,ilz), 0)
endif
else if (proc0) then
call receive (ltmp, proc_id(lz_lo, ilz))
call receive (ftmp, proc_id(lz_lo, ilz))
endif
call barrier
do il=1,te
etmp = cabs(ltmp(il) - ftmp(il))
if (etmp > emax) then
emax = etmp
ltmax = cabs(ltmp(il))
ikmax = ik
itmax = it
iemax = ie
ismax = is
ilmax = il
igmax = ig
end if
etot = etot + etmp
edenom = edenom + cabs(ltmp(il))
end do
end do
eavg = etot/edenom
emax = emax/ltmax
if (proc0) then
write(unit,"((1x,e13.6),6(i8),2(1x,e13.6))") time, &
igmax, ikmax, itmax, iemax, ilmax, ismax, emax, eavg
if (last) then
call close_output_file (unit)
end if
end if
deallocate (lcoll, fdcoll, glze, ltmp, ftmp)
end subroutine collision_error