Calculate the time-averaged complex frequency, check convergence criterion, and numerical instability criterion.
Time average is over previous navg timesteps
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | istep |
The current timestep |
||
| logical, | intent(inout) | :: | exit |
Should the simulation exit?
FIXME: This could be |
||
| complex, | intent(out), | dimension (:,:) | :: | omegaavg |
The time-averaged complex frequency |
|
| logical, | intent(in), | optional | :: | debopt |
If true, write some debug messages |
subroutine get_omegaavg (istep, exit, omegaavg, debopt)
use kt_grids, only: ntheta0, naky
use gs2_time, only: code_dt
implicit none
!> The current timestep
integer, intent (in) :: istep
!> Should the simulation exit?
!> FIXME: This could be `intent(out)` only
logical, intent (in out) :: exit
!> The time-averaged complex frequency
complex, dimension (:,:), intent (out) :: omegaavg
!> If true, write some debug messages
logical, intent (in), optional :: debopt
logical :: debug=.false.
if (debug) write(6,*) "get_omeaavg: allocate domega"
if (.not. allocated(domega)) allocate(domega(navg,ntheta0,naky))
if (present(debopt)) debug=debopt
if (debug) write(6,*) "get_omeaavg: start"
if (debug) write(6,*) "get_omeaavg: set omegahist"
! Get and store the current instantaneous omega
omegahist(mod(istep, navg), :, :) = calculate_instantaneous_omega(igomega, omegatol)
if (debug) write(6,*) "get_omeaavg: set omegahist at istep = 0"
!During initialisation fieldnew==field but floating error can lead to finite omegahist
!Force omegahist=0 here to avoid erroneous values.
!Could think about forcing omegahist=0 where abs(omegahist)<tol
if(istep.eq.0) omegahist(:,:,:)=0.0
if (debug) write(6,*) "get_omeaavg: set omegaavg"
omegaavg = sum(omegahist/real(navg),dim=1)
if (debug) write(6,*) "get_omegaavg: omegaavg=",omegaavg
if (istep > navg) then
domega = spread(omegaavg,1,navg) - omegahist
if (all(sqrt(sum(abs(domega)**2/real(navg),dim=1)) &
.le. min(abs(omegaavg),1.0)*omegatol)) &
then
if (write_ascii) write (out_unit, "('*** omega converged')")
exit = .true. .and. exit_when_converged
end if
if (any(abs(omegaavg)*code_dt > omegatinst)) then
if (write_ascii) write (out_unit, "('*** numerical instability detected')")
exit = .true.
end if
end if
if (debug) write(6,*) "get_omegaavg: done"
end subroutine get_omegaavg