fm_get_condition_numbers Subroutine

private subroutine fm_get_condition_numbers(self)

Fetch all condition numbers to proc0 and report on worrying values as required.

Type Bound

fieldmat_type

Arguments

Type IntentOptional Attributes Name
class(fieldmat_type), intent(in) :: self

Contents


Source Code

  subroutine fm_get_condition_numbers(self)
    use file_utils, only: error_unit
    use gs2_time, only: code_dt
    use mp, only: proc0, sum_reduce
    use unit_tests, only: debug_message
    implicit none
    type ky_condition_type
       real, dimension(:), allocatable ::  conditions
    end type ky_condition_type
    class(fieldmat_type), intent(in) :: self
    integer, parameter :: verb = 1
    type(ky_condition_type), dimension(:), allocatable :: ky_conditions
    integer :: ik, is
    character(len = 100) :: condition_number_message
    real :: condition_number
    allocate(ky_conditions(self%naky))

    ! Note that here we communciate over comm_world. It is likely that we
    ! could instead reduce individual `ky_conditions(:)%conditions` arrays
    ! amongst just the supercell heads of the ky block and then reduce
    ! across the ky heads to proc0. Currently we don't worry about this as
    ! we just do this once per initialisation, but this could become expensive
    ! when running at large scale?

    do ik = 1, self%naky
       allocate(ky_conditions(ik)%conditions(self%kyb(ik)%nsupercell))
       ky_conditions(ik)%conditions = 0.0
       do is = 1, self%kyb(ik)%nsupercell
          if (self%kyb(ik)%supercells(is)%is_head) then
             ky_conditions(ik)%conditions(is) = self%kyb(ik)%supercells(is)%condition_number
          end if
       end do
       call sum_reduce(ky_conditions(ik)%conditions, 0)
    end do

    if (proc0) then
       do ik = 1, self%naky
          do is = 1, self%kyb(ik)%nsupercell
             condition_number = ky_conditions(ik)%conditions(is)
             write(condition_number_message,'("Condition number for response with ik = ",I0,", is = ",I0," and dt = ",E14.6E3," is ",E14.6E3)') ik, is, code_dt, condition_number
             call debug_message(verb, "fields_local::fm_get_condition_numbers "//condition_number_message)
             if (condition_number > condition_number_warning_limit) then
                write(error_unit(), '("Warning: ",A)') trim(condition_number_message)
             end if

             ! In debug mode always display the condition number message to screen
             if (debug) then
                write(*,'(A)') trim(condition_number_message)
             end if
          end do
       end do
    end if

    do ik = 1, self%naky
       deallocate(ky_conditions(ik)%conditions)
    end do
    deallocate(ky_conditions)
  end subroutine fm_get_condition_numbers