Calculate and write the various diagnostic quantities.
This is the main routine for the "old" diagnostics.
The rest of the routine only happens every nwrite steps
The following large section could do with being moved to separate routines but at the moment its all quite interlinked which makes this hard.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | istep |
Current timestep |
||
| logical, | intent(out) | :: | exit |
If true, simulation should stop (for example, because it has converged) |
||
| logical, | intent(in), | optional | :: | debopt |
If true, turn on some debug prints |
|
| logical, | intent(in), | optional | :: | force |
subroutine loop_diagnostics (istep, exit, debopt, force)
use, intrinsic :: iso_fortran_env, only: output_unit
use species, only: nspec, spec, has_electron_species
use kt_grids, only: naky, ntheta0
use run_parameters, only: nstep, has_phi, has_apar, has_bpar
use fields_arrays, only: phinew, bparnew
use collisions, only: ncheck, vary_vnew
use mp, only: proc0, broadcast
use gs2_time, only: user_time
use gs2_io, only: nc_qflux, nc_vflux, nc_pflux, nc_pflux_tormom, nc_exchange, nc_final_fields, nc_sync, get_netcdf_movie_file_id
use nonlinear_terms, only: nonlin
use antenna, only: antenna_w
use optionals, only: get_option_with_default
use diagnostics_printout, only: write_phi_fluxes_to_unit, write_fluxes_to_unit
use diagnostics_fluxes, only: common_calculate_fluxes, qheat, qmheat, qbheat, &
pflux, vflux, vflux_par, vflux_perp, pflux_tormom, vflux0, vflux1, pmflux, vmflux, &
pbflux, vbflux, exchange
use diagnostics_velocity_space, only: collision_error
use optionals, only: get_option_with_default
use warning_helpers, only: is_not_zero, not_exactly_equal
implicit none
!> Current timestep
integer, intent (in) :: istep
!> If true, simulation should stop (for example, because it has converged)
logical, intent (out) :: exit
!> If true, turn on some debug prints
logical, intent (in), optional:: debopt
logical, intent(in), optional :: force
real, dimension (ntheta0, naky) :: phitot
real :: phi2, apar2, bpar2
real :: t
integer :: ik, it, is, write_mod
complex, save :: wtmp_new !This shouldn't need to be given the save property
real, dimension (ntheta0, nspec) :: x_qmflux
real, dimension (nspec) :: heat_fluxes, part_fluxes, mom_fluxes, parmom_fluxes, perpmom_fluxes, part_tormom_fluxes
real, dimension (nspec) :: mheat_fluxes, mpart_fluxes, mmom_fluxes
real, dimension (nspec) :: bheat_fluxes, bpart_fluxes, bmom_fluxes
real, dimension (nspec) :: energy_exchange
real, dimension (nspec) :: heat_par, heat_perp
real, dimension (nspec) :: mheat_par, mheat_perp
real, dimension (nspec) :: bheat_par, bheat_perp
real :: hflux_tot, zflux_tot, vflux_tot
logical :: do_force, debug
real, save :: t_old = 0.
integer, save :: istep_last = -1
debug = get_option_with_default(debopt, .false.)
!Set the current time
t = user_time
!Always skip if we've already done this step
if (istep == istep_last) return
do_force = get_option_with_default(force, .false.)
exit = .false.
call do_get_omega(istep,debug,exit)
if (write_heating) call heating (istep, navg, h, hk)
if (make_movie .and. mod(istep,nmovie) == 0) call do_write_movie(t)
if (write_gyx .and. mod(istep,nmovie) == 0) call do_write_fyx (yxdist_unit, phinew, bparnew)
if (vary_vnew) then
write_mod = mod(istep,ncheck)
else
write_mod = mod(istep,nwrite)
end if
if (write_verr .and. write_mod == 0) call do_write_verr
if (debug) write(6,*) "loop_diagnostics: call update_time"
!########################################################
!The rest of the routine only happens every nwrite steps
!########################################################
if (mod(istep,nwrite) /= 0 .and. .not. exit .and. .not. do_force) return
! If we get here then we're doing a full set of diagnostics
! so store the step
istep_last = istep
if (write_g) call do_write_f (dist_unit)
! Note this also returns phi2, apar2, bpar2 and phitot for other diagnostics
if (proc0) call do_write_ncloop(t, istep, phi2, apar2, bpar2, phitot)
if (print_line) call do_print_line(phitot)
if (write_moments) call do_write_moments
if (write_cross_phase .and. has_electron_species(spec) .and. write_ascii) call do_write_crossphase(t)
!###########################
! The following large section could do with being moved to separate routines but
! at the moment its all quite interlinked which makes this hard.
!###########################
!Zero various arrays which may or may not get filled
part_fluxes = 0.0 ; mpart_fluxes = 0.0 ; bpart_fluxes = 0.0
heat_fluxes = 0.0 ; mheat_fluxes = 0.0 ; bheat_fluxes = 0.0
mom_fluxes = 0.0 ; mmom_fluxes = 0.0 ; bmom_fluxes = 0.0
part_tormom_fluxes = 0.0
energy_exchange = 0.0
if (debug) write(6,*) "loop_diagnostics: -1"
if (write_any_fluxes) then
call common_calculate_fluxes()
if (proc0) then
if (has_phi) then
do is = 1, nspec
call get_volume_average (qheat(:,:,is,1), heat_fluxes(is))
call get_volume_average (qheat(:,:,is,2), heat_par(is))
call get_volume_average (qheat(:,:,is,3), heat_perp(is))
call get_volume_average (pflux(:,:,is), part_fluxes(is))
call get_volume_average (pflux_tormom(:,:,is), part_tormom_fluxes(is))
call get_volume_average (vflux(:,:,is), mom_fluxes(is))
call get_volume_average (vflux_par(:,:,is), parmom_fluxes(is))
call get_volume_average (vflux_perp(:,:,is), perpmom_fluxes(is))
call get_volume_average (exchange(:,:,is), energy_exchange(is))
end do
end if
if (has_apar) then
do is = 1, nspec
call get_volume_average (qmheat(:,:,is,1), mheat_fluxes(is))
call get_surf_average (qmheat(:,:,is,1), x_qmflux(:,is))
call get_volume_average (qmheat(:,:,is,2), mheat_par(is))
call get_volume_average (qmheat(:,:,is,3), mheat_perp(is))
call get_volume_average (pmflux(:,:,is), mpart_fluxes(is))
call get_volume_average (vmflux(:,:,is), mmom_fluxes(is))
end do
end if
if (has_bpar) then
do is = 1, nspec
call get_volume_average (qbheat(:,:,is,1), bheat_fluxes(is))
call get_volume_average (qbheat(:,:,is,2), bheat_par(is))
call get_volume_average (qbheat(:,:,is,3), bheat_perp(is))
call get_volume_average (pbflux(:,:,is), bpart_fluxes(is))
call get_volume_average (vbflux(:,:,is), bmom_fluxes(is))
end do
end if
pflux_avg = pflux_avg + (part_fluxes + mpart_fluxes + bpart_fluxes)*(t-t_old)
qflux_avg = qflux_avg + (heat_fluxes + mheat_fluxes + bheat_fluxes)*(t-t_old)
vflux_avg = vflux_avg + (mom_fluxes + mmom_fluxes + bmom_fluxes)*(t-t_old)
if (write_heating) heat_avg = heat_avg + h%imp_colls*(t-t_old)
! t_old = t
end if
!These are for trinity -- can we detect if this is a trinity run
!and avoid these broadcasts when not trinity?
call broadcast (pflux_avg)
call broadcast (qflux_avg)
call broadcast (vflux_avg)
if (write_heating) call broadcast (heat_avg)
end if
if (proc0 .and. print_flux_line) then
if (has_phi) then
call write_phi_fluxes_to_unit(output_unit, t, phi2, heat_fluxes, energy_exchange)
end if
if (has_apar) then
call write_fluxes_to_unit(output_unit, t, apar2, "apar", mheat_fluxes)
end if
if (has_bpar) then
call write_fluxes_to_unit(output_unit, t, bpar2, "bpar", bheat_fluxes)
end if
end if
! Check for convergence
if (nonlin .and. use_nonlin_convergence) call check_nonlin_convergence(istep, heat_fluxes(1), exit)
if (debug) write(6,*) "loop_diagnostics: -1"
if (.not. (write_any .or. dump_any)) then
! We have to increment the number of outputs written to file
! before leaving. Usually we do this at the end of the routine
! so we must make sure we do this here before leaving early.
nout = nout + 1
return
end if
if (debug) write(output_unit, *) "loop_diagnostics: -2"
if (proc0 .and. write_any) then
if (write_ascii) then
write (unit=out_unit, fmt=*) 'time=', t
if (write_heating) call do_write_heating(t, heat_unit, heat_unit2, h)
if (write_jext) call do_write_jext(t,istep)
end if
hflux_tot = 0.
zflux_tot = 0.
vflux_tot = 0.
if (has_phi) then
hflux_tot = sum(heat_fluxes)
vflux_tot = sum(mom_fluxes)
zflux_tot = sum(part_fluxes*spec%z)
end if
if (has_apar) then
hflux_tot = hflux_tot + sum(mheat_fluxes)
vflux_tot = vflux_tot + sum(mmom_fluxes)
zflux_tot = zflux_tot + sum(mpart_fluxes*spec%z)
end if
if (has_bpar) then
hflux_tot = hflux_tot + sum(bheat_fluxes)
vflux_tot = vflux_tot + sum(bmom_fluxes)
zflux_tot = zflux_tot + sum(bpart_fluxes*spec%z)
end if
if (write_flux_line .and. write_ascii) then
if (has_phi) then
call write_phi_fluxes_to_unit(out_unit, t, phi2, &
heat_fluxes, energy_exchange, part_fluxes, mom_fluxes)
end if
if (write_lorentzian) then
wtmp_new = antenna_w()
if (is_not_zero(real(wtmp_old)) .and. not_exactly_equal(wtmp_new, wtmp_old)) &
write (unit=out_unit, fmt="('w= ',e17.10, ' amp= ',e17.10)") real(wtmp_new), sqrt(2.*apar2)
wtmp_old = wtmp_new
end if
if (has_apar) then
call write_fluxes_to_unit(out_unit, t, apar2, "apar", mheat_fluxes, mpart_fluxes)
end if
if (has_bpar) then
call write_fluxes_to_unit(out_unit, t, bpar2, "bpar", bheat_fluxes, bpart_fluxes)
end if
write (unit=out_unit, fmt="('t= ',e17.10,' h_tot= ',e11.4, ' z_tot= ',e11.4)") t, hflux_tot, zflux_tot
end if
if (write_fluxes) then
call nc_qflux (get_netcdf_file_id(), nout, qheat(:,:,:,1), qmheat(:,:,:,1), qbheat(:,:,:,1), &
heat_par, mheat_par, bheat_par, &
heat_perp, mheat_perp, bheat_perp, &
heat_fluxes, mheat_fluxes, bheat_fluxes, x_qmflux, hflux_tot)
call nc_exchange (get_netcdf_file_id(), nout, exchange, energy_exchange)
call nc_vflux (get_netcdf_file_id(), nout, vflux, vmflux, vbflux, &
mom_fluxes, mmom_fluxes, bmom_fluxes, vflux_tot, &
vflux_par, vflux_perp, vflux0, vflux1)
call nc_pflux (get_netcdf_file_id(), nout, pflux, pmflux, pbflux, &
part_fluxes, mpart_fluxes, bpart_fluxes, zflux_tot)
end if
if (write_ascii) then
do ik = 1, naky
do it = 1, ntheta0
if (write_line) call do_write_line(t,it,ik,phitot(it,ik))
if (write_omega) call do_write_omega(it,ik)
if (write_omavg) call do_write_omavg(it,ik)
end do
enddo
end if
endif
if (write_cerr) call collision_error(cres_unit, phinew, bparnew)
if (write_nl_flux_dist) call do_write_nl_flux_dist(get_netcdf_file_id(), nout)
if (write_kinetic_energy_transfer) call do_write_kinetic_energy_transfer(get_netcdf_file_id(), nout)
if (write_symmetry) call do_write_symmetry(get_netcdf_file_id(), nout)
if (write_correlation) call do_write_correlation(get_netcdf_file_id(), nout)
if (write_correlation_extend .and. istep > nstep/4 .and. mod(istep,nwrite_mult*nwrite)==0) &
call do_write_correlation_extend(get_netcdf_file_id(), t, t_old)
if (write_parity) call do_write_parity(t, parity_unit, write_ascii)
if (write_avg_moments) call do_write_avg_moments
if (write_full_moments_notgc) call do_write_full_moments_notgc
if (dump_check1) call do_write_dump_1(t)
if (dump_check2) call do_write_dump_2(t)
if (dump_fields_periodically .and. mod(istep,10*nwrite) == 0) call do_dump_fields_periodically(t)
!Now sync the data to file (note doesn't actually sync every call)
if (proc0) call nc_sync(get_netcdf_file_id(), nout, get_netcdf_movie_file_id(), nout_movie, nc_sync_freq)
!Increment loop counter
nout = nout + 1
if (write_ascii .and. mod(nout, 10) == 0 .and. proc0) call flush_files
!Update time
t_old = t
if (debug) write(6,*) "loop_diagnostics: done"
end subroutine loop_diagnostics