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 build_config, only: gs2_has_lowflow
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 dist_fn, only: collision_error, write_f, write_fyx
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 parameter_scan_arrays, only: scan_hflux => hflux_tot, scan_momflux => momflux_tot
use parameter_scan_arrays, only: scan_phi2_tot => phi2_tot, scan_nout => nout
use parameter_scan, only: scanning
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
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) :: lfmom_fluxes, vflux1_avg ! low-flow correction to turbulent momentum 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
real, save :: t_old = 0.
logical :: last = .false.
logical:: debug=.false.
logical :: do_force
integer, save :: istep_last = -1
if (present(debopt)) debug = debopt
!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, h, hk)
if (make_movie .and. mod(istep,nmovie) == 0) call do_write_movie(t)
if (write_gyx .and. mod(istep,nmovie) == 0) call write_fyx (phinew,bparnew,last)
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 write_f (last)
! 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))
#ifdef LOWFLOW
call get_volume_average (vflux0(:,:,is), lfmom_fluxes(is))
call get_volume_average (vflux1(:,:,is), vflux1_avg(is))
#endif
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, &
lfmom_fluxes, vflux1_avg, gs2_has_lowflow)
end if
if (write_lorentzian) then
wtmp_new = antenna_w()
if (real(wtmp_old) /= 0. .and. 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)
! Update the target array in parameter_scan_arrays
!Do you really want the scan parameter to rely on a diagnostic flag?
scan_hflux(mod(nout-1,nstep/nwrite+1)+1) = hflux_tot
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)
! Update the target array in parameter_scan_arrays
!Do you really want the scan parameter to rely on a diagnostic flag?
scan_momflux(mod(nout-1,nstep/nwrite+1)+1) = vflux_tot
call nc_pflux (get_netcdf_file_id(), nout, pflux, pmflux, pbflux, &
part_fluxes, mpart_fluxes, bpart_fluxes, zflux_tot)
end if
if (write_pflux_tormom) call nc_pflux_tormom (get_netcdf_file_id(), nout, pflux_tormom, part_tormom_fluxes)
scan_phi2_tot(mod(nout-1,nstep/nwrite+1)+1) = phi2
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 (scanning) call bcast_scan_parameter(scan_hflux, scan_momflux, scan_phi2_tot)
if (write_cerr) call collision_error(phinew, bparnew, last)
if (write_nl_flux_dist) call do_write_nl_flux_dist(get_netcdf_file_id(), nout)
! In new diagnostics, these two are both called for `write_symmetry`
if (write_symmetry) call do_write_symmetry(get_netcdf_file_id(), nout)
if (write_pflux_sym) call do_write_pflux_sym(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)
! Update the counter in parameter_scan_arrays
scan_nout = nout
!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