init_fields_matrixlocal Subroutine

private subroutine init_fields_matrixlocal(tuning_in)

Routine use to initialise field matrix data !!!!!!!! NOTE: All of the above should be a one off setup cost unless something fundamental changes with the parallel layout (a change in nproc/layout etc.) The remaining two tasks (populate+prepare) must be repeated each time the relevant physics and numerical parameters change. In practice this just means when the time step changes. !!!!!!!!

Arguments

Type IntentOptional Attributes Name
logical, optional :: tuning_in

Contents


Source Code

  subroutine init_fields_matrixlocal(tuning_in)
    use mp, only: proc0, mp_abort, land_allreduce
    use file_utils, only: error_unit
    use optionals, only: get_option_with_default
    implicit none
    logical, optional :: tuning_in
    logical :: tuning
    logical :: response_was_read
    real :: ts,te

    !Exit if we've already initialised
    if(initialised) return

    tuning = get_option_with_default(tuning_in, .false.)

    if(.not.reinit)then
       !Use the fieldmat init routine to create fieldmat
       !structure and fill some of the basic variables
       !if(debug)call barrier
       if(proc0.and.debug) then
          write(dlun,'("Initialising fieldmat.")')
          call cpu_time(ts)
       endif
       call fieldmat%init
       !if(debug)call barrier
       if(proc0.and.debug) then 
          call cpu_time(te)
          write(dlun,'("--Done in ",F12.4," units")') te-ts
       endif

       !Now initialise the parallel decomposition object
       !Note: We don't actually decompose until we make some
       !of the primary subcom so we know how many procs are 
       !available to each object etc.
       !if(debug)call barrier
       if(proc0.and.debug) then
          write(dlun,'("Initialising pc.")')
          call cpu_time(ts)
       endif
       call pc%init
       !if(debug)call barrier
       if(proc0.and.debug) then 
          call cpu_time(te)
          write(dlun,'("--Done in ",F12.4," units")') te-ts
       endif

       !Now use the pc object to set field object locality
       !and make the primary subcommunicators
       !if(debug)call barrier
       if(proc0.and.debug) then
          write(dlun,'("Setting initial locality and making primary subcommunicators.")')
          call cpu_time(ts)
       endif
       call fieldmat%set_is_local
       call fieldmat%make_subcom_1
       !if(debug)call barrier
       if(proc0.and.debug) then 
          call cpu_time(te)
          write(dlun,'("--Done in ",F12.4," units")') te-ts
       endif

       !Now we can setup the decomposition, which basically
       !consists of setting the row limit of row block objects
       !if(debug)call barrier
       if(proc0.and.debug) then
          write(dlun,'("Initialising decomposition.")')
          call cpu_time(ts)
       endif
       call pc%make_decomp
       !if(debug)call barrier
       if(proc0.and.debug) then 
          call cpu_time(te)
          write(dlun,'("--Done in ",F12.4," units")') te-ts
       endif

       !Now set fieldmat locality and allocate data arrays
       !if(debug)call barrier
       if(proc0.and.debug) then
          write(dlun,'("Setting locality and allocating.")')
          call cpu_time(ts)
       endif
       call fieldmat%set_is_local
       call fieldmat%allocate
       !if(debug)call barrier
       if(proc0.and.debug) then 
          call cpu_time(te)
          write(dlun,'("--Done in ",F12.4," units")') te-ts
       endif

       !Need to make cell/supercell level sub-communicators
       !at some point before we try to invert.
       !Can make them at any point after locality defined

       !/Note that to split a sub-comm everybody in the original
       !communicator needs to make the call. We can either split
       !all groups off from WORLD or we can take "nested subsets"
       !Note sure if there is any perfomance benefit in the actual
       !comms (i'd guess not) but it should be faster to split smaller
       !groups and can do it in parallel so the nested subset approach
       !may be better.
       !if(debug)call barrier
       if(proc0.and.debug)then
          write(dlun,'("Creating the secondary subcoms.")')
          call cpu_time(ts)
       endif
       call fieldmat%make_subcom_2

       call fieldmat%calc_sc_heads

       !if(debug)call barrier
       if(proc0.and.debug) then 
          call cpu_time(te)
          write(dlun,'("--Done in ",F12.4," units")') te-ts
       endif
    endif
!!!!!!!!!!
!!NOTE: All of the above should be a one off setup cost unless something
!!fundamental changes with the parallel layout (a change in nproc/layout etc.)
!!The remaining two tasks (populate+prepare) must be repeated each time the
!!relevant physics and numerical parameters change. In practice this just means
!!when the time step changes.
!!!!!!!!!!
    if(.not. tuning) then

       response_was_read = .false.

       !Now fill the fieldmat with data
       if(read_response)then
          if(proc0.and.debug)then
             write(dlun,'("Populating from dump.")')
             call cpu_time(ts)
          endif
          call read_response_from_file_local(could_read = response_was_read)
          if(proc0.and.debug) then
             call cpu_time(te)
             write(dlun,'("--Done in ",F12.4," units")') te-ts
          endif

          ! Ensure we reduce response_was_read across all processors so if _any_
          ! processors couldn't read the matrices we recalculate on all processors.
          call land_allreduce(response_was_read)
          if((.not. response_was_read) .and. proc0) write(error_unit(), &
               '("Could not find response file so reverting to calculation.")')

       end if

       !If we weren't able to read the response file, or didn't try to,
       !then we need to calculate the response matrix instead
       if (.not. response_was_read) then
          !if(debug)call barrier
          if(proc0.and.debug)then
             write(dlun,'("Populating.")')
             call cpu_time(ts)
          endif
          call fieldmat%populate
          !if(debug)call barrier
          if(proc0.and.debug) then
             call cpu_time(te)
             write(dlun,'("--Done in ",F12.4," units")') te-ts
          endif
          
          !Now prepare the response matrices
          !Note: Currently prepare just means invert
          !but at some point we may wish to LU decompose
          !or other approaches etc.
          !if(debug)call barrier
          if(proc0.and.debug) then
             write(dlun,'("Preparing.")')
             call cpu_time(ts)
          endif
          call fieldmat%prepare
          !if(debug)call barrier
          if(proc0.and.debug) then 
             call cpu_time(te)
             write(dlun,'("--Done in ",F12.4," units")') te-ts
          endif
          call fieldmat%get_condition_numbers
       endif

       !Dump response to file
       if(dump_response)then
          if(proc0.and.debug)then
             write(dlun,'("Dumping to file.")')
             call cpu_time(ts)
          endif
          call dump_response_to_file_local
          if(proc0.and.debug) then
             call cpu_time(te)
             write(dlun,'("--Done in ",F12.4," units")') te-ts
          endif
       endif

    endif
    !Now write debug data
!    call fieldmat%write_debug_data
  end subroutine init_fields_matrixlocal