minimal_model_v5.f90

Go to the documentation of this file.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    minimal_model_comms just sets up communicators for empire for
!    testing purposes
!
!The MIT License (MIT)
!
!Copyright (c) 2015 Philip A. Browne
!
!Permission is hereby granted, free of charge, to any person obtaining a copy
!of this software and associated documentation files (the "Software"), to deal
!in the Software without restriction, including without limitation the rights
!to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
!copies of the Software, and to permit persons to whom the Software is
!furnished to do so, subject to the following conditions:
!
!The above copyright notice and this permission notice shall be included in all
!copies or substantial portions of the Software.
!
!THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
!AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
!OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
!SOFTWARE.
!
!Email: p.browne@reading.ac.uk
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program minimal_model_v5
  implicit none

  include 'mpif.h'
  integer :: mpi_err
  integer :: world_rank
  integer :: world_size
  integer, parameter :: mdl_num_proc=16
  integer, parameter :: ensemble_members_per_proc=2
  integer :: mdl_mpi_comm
  integer :: temp_mdls_comm
  integer :: temp_mdls_rank
  integer :: mdl_rank
  integer :: da
  integer, parameter :: rk = kind(1.0d0)
  real(kind=rk), allocatable, dimension(:) :: state_vector
  integer :: state_dim
  integer :: i,j
  integer :: cpl_root
  integer :: instance_rank
  integer :: particle_rank
  integer :: nda
  integer :: nens
  integer :: temp_mdls_size
  integer :: temp_cpl_comm
  integer :: first_ptcl
  integer :: final_ptcl
!  integer :: cpl_mpi_comm
  integer, dimension(ensemble_members_per_proc) :: cpl_mpi_comms
  integer :: null_mpi_comm
  real(kind=rk), dimension(0) :: send_null
  integer :: total_timesteps
  integer :: tag

  call mpi_init(mpi_err)
      if(mpi_err .eq. 0) then    
       print*,'mpi_init successful'
    else
       print*,'mpi_init unsuccessful'
    end if


  call mpi_comm_rank(mpi_comm_world,world_rank,mpi_err)
    if(mpi_err .eq. 0) then    
       print*,'mpi_comm_rank successful'
       print*,'world_rank = ',world_rank
    else
       print*,'mpi_comm_rank unsuccessful'
    end if


  call mpi_comm_size(mpi_comm_world,world_size,mpi_err)
      if(mpi_err .eq. 0) then    
       print*,'mpi_comm_size successful'
       print*,'world_size = ',world_size
    else
       print*,'mpi_comm_size unsuccessful'
    end if


  cpl_root = world_size-1
  print*,'rank = ',world_rank,' on mpi_comm_world which has size ',world_size
  da = 0


  call mpi_allreduce(mdl_num_proc,i,1,mpi_integer,mpi_max&
       &,mpi_comm_world,mpi_err)
  if(mpi_err .eq. 0) then    
       print*,'mpi_allreduce successful'
       print*,'i = ',i
    else
       print*,'mpi_allreduce unsuccessful'
    end if


  call mpi_allreduce(ensemble_members_per_proc,i,1,mpi_integer,mpi_max&
       &,mpi_comm_world,mpi_err)
  if(mpi_err .eq. 0) then    
     print*,'mpi_allreduce successful'
     print*,'i = ',i
  else
     print*,'mpi_allreduce unsuccessful'
  end if


  print*,'first split happening with the inputs: ', mpi_comm_world,da,world_rank,temp_mdls_comm&
       &,mpi_err
  call mpi_comm_split(mpi_comm_world,da,world_rank,temp_mdls_comm&
       &,mpi_err)
    print*,'first split happened with the outputs: ', mpi_comm_world&
         &,da,world_rank,temp_mdls_comm,mpi_err
      if(mpi_err .eq. 0) then    
       print*,'mpi_comm_split successful: temp_mdls_comm created: ',temp_mdls_comm
    else
       print*,'mpi_comm_split unsuccessful: temp_mdls_comm not created'
    end if
  
    call mpi_comm_size(temp_mdls_comm,temp_mdls_size,mpi_err)
    if(mpi_err .eq. 0) then    
     print*,'mpi_comm_size successful'
     print*,'temp_mdls_size = ',temp_mdls_size
  else
     print*,'mpi_comm_size unsuccessful'
  end if
  

  if(mod(temp_mdls_size,mdl_num_proc) .ne. 0) then
     print*,'MINIMAL MODEL LAUNCH ERROR.'
     print*,'MUST LAUNCH A MULTIPLE OF ',mdl_num_proc,' copies of the &
          &model'
     stop
  end if


  nda = world_size-temp_mdls_size
  print*,'nda = ',nda
  if(nda .lt. 1) then
     print*,'MINIMAL MODEL COMMS v5 ERROR: nda is less than 1.'
     print*,'Make sure you launch with a DA CODE'
     stop
  end if



  nens = (temp_mdls_size/mdl_num_proc)*ensemble_members_per_proc
  print*,'nens = ',nens
  call mpi_comm_rank(temp_mdls_comm,temp_mdls_rank,mpi_err)
    if(mpi_err .eq. 0) then    
       print*,'mpi_comm_rank successful'
       print*,'temp_mdls_rank = ',temp_mdls_rank
    else
       print*,'mpi_comm_rank unsuccessful'
    end if


  instance_rank = temp_mdls_rank/mdl_num_proc

  call mpi_comm_split(temp_mdls_comm,instance_rank,temp_mdls_rank&
       &,mdl_mpi_comm,mpi_err)
  if(mpi_err .eq. 0) then    
       print*,'mpi_comm_split successful: mdl_mpi_comm created ',mdl_mpi_comm
    else
       print*,'mpi_comm_split unsuccessful: mdl_mpi_comm not created'
    end if
    


  call mpi_comm_rank(mdl_mpi_comm,mdl_rank,mpi_err)
    if(mpi_err .eq. 0) then    
       print*,'mpi_comm_rank successful'
       print*,'mdl_rank = ',mdl_rank
    else
       print*,'mpi_comm_rank unsuccessful'
    end if


  cpl_root = nda*instance_rank/(nens/ensemble_members_per_proc)
  print*,'cpl_root = ',cpl_root

  if(cpl_root .lt. 0) then
     print*,'MINIMAL MODEL COMMS v5 ERROR: cpl_root is less than 0.'
     print*,'Make sure you launch with a DA CODE'
     stop
  end if


  first_ptcl = -((-cpl_root)*nens/nda)
  final_ptcl = -((-cpl_root-1)*nens/nda)-1

  first_ptcl = ceiling(real(cpl_root)*real(nens)/real(nda))
  final_ptcl = ceiling(real(cpl_root+1)*real(nens)/real(nda))-1


  print*,'range of particles = ',first_ptcl,final_ptcl

 
  !check to see if mod(n_mdl_instances,npfs) == 0:
  if(mod(temp_mdls_size/mdl_num_proc,nda) .ne. 0 ) then
     print*,'EMPIRE MESSAGE: SEQUENTIAL SPLITTING OF MPI_COMM_WORLD'
     j = 0
     
     do i = 0,nens-1
        print*,'i = ',i,' particle_rank = ',particle_rank
        print*,'instance_rank*ensemble_members_per_proc = ',instance_rank*ensemble_members_per_proc
        print*,'(instance_rank+1)*ensemble_members_per_proc -1 = '&
             &,(instance_rank+1)*ensemble_members_per_proc -1
        if(i .ge. instance_rank*ensemble_members_per_proc .and. &
             i .le. (instance_rank+1)*ensemble_members_per_proc -1) then
           j = j + 1
           call mpi_comm_split(mpi_comm_world,1,temp_mdls_rank&
                &,cpl_mpi_comms(j),mpi_err)
           print*,'created cpl_mpi_comms(j) with number',cpl_mpi_comms(j)
        else
           print*,'doing null splitting'
           call mpi_comm_split(mpi_comm_world,0,temp_mdls_rank&
                &,null_mpi_comm,mpi_err)
           print*,'created mpi_comm_null'
           call mpi_comm_free(null_mpi_comm,mpi_err)
           print*,'freed up null_mpi_comm'
        end if


     end do
     

     
  else
        print*,''
       print*,''
       print*,''
       print*,''
       print*,''

       print*,'doing first split based on pfrank',cpl_root
       print*,'doinf first split...key = ',temp_mdls_rank
       temp_cpl_comm = 0
       print*,mpi_comm_world,cpl_root,temp_mdls_rank&
          &,temp_cpl_comm,mpi_err
     call mpi_comm_split(mpi_comm_world,cpl_root,temp_mdls_rank&
          &,temp_cpl_comm,mpi_err)
     if(mpi_err .eq. 0) then    
        print*,'mpi_comm_split successful: temp_cpl_comm created ',temp_cpl_comm
     else
        print*,'mpi_comm_split unsuccessful: temp_cpl_comm not created'
     end if
     
     print*,'finished first split'
     print*,''
       print*,''
       print*,''
       print*,''
       print*,''

     
     j = 0

  do i = first_ptcl,final_ptcl
     print*,'i = ',i,' particle_rank = ',particle_rank
        if(i .ge. instance_rank*ensemble_members_per_proc .and. &
             i .le. (instance_rank+1)*ensemble_members_per_proc -1) then
           j = j + 1
        call mpi_comm_split(temp_cpl_comm,1,temp_mdls_rank&
                &,cpl_mpi_comms(j),mpi_err)
           print*,'created cpl_mpi_comms(j) with number',cpl_mpi_comms(j)
     else
        print*,'doing null splitting'
        call mpi_comm_split(temp_cpl_comm,0,temp_mdls_rank&
             &,null_mpi_comm,mpi_err)
        print*,'created mpi_comm_null'
        call mpi_comm_free(null_mpi_comm,mpi_err)
        print*,'freed up null_mpi_comm'
     end if


  end do




  end if


  cpl_root = mdl_num_proc



  select case(mdl_rank)
  case(0)
     state_dim = 1
  case(1)
     state_dim = 3
  case(2)
     state_dim = 2
  case(3)
     state_dim = 5
  case(4)
     state_dim = 1
  case(5,6,7,8,9,10,11,12,13,14,15)
     state_dim = 2
  case default
     print*,'it was at this point, model realised, he fucked up'
     stop
  end select
  allocate(state_vector(state_dim))




  state_vector = 10*mdl_rank + (/ (real(i,rk), i = 1,state_dim) /)

  print*,'state_vector = '
  print*,state_vector

  do j = 1,ensemble_members_per_proc
     print*,'doing a gather on cpl_mpi_comm(j)',cpl_mpi_comms(j),cpl_root
  call mpi_gather(state_dim,1,mpi_integer,state_dim&
       &,1,mpi_integer,cpl_root,cpl_mpi_comms(j),mpi_err)
  print*,'finished the gather on cpl_mpi_comms(j)'
  end do





















  print*,'Reading total_timesteps from file timesteps: '
  open(11,file='timesteps',action='read',status='old')
  read(11,*) total_timesteps
  close(11)

  do j = 1,ensemble_members_per_proc
     !send the state to da code with mpi_gatherv
     call mpi_gatherv(state_vector,state_dim,mpi_double_precision,state_vector&
          &,state_dim,state_dim,mpi_double_precision,cpl_root&
          &,cpl_mpi_comms(j),mpi_err)
  end do

  do j = 1,ensemble_members_per_proc
     !get the state back from da code with mpi_gatherv
     call mpi_scatterv(send_null,0,0,mpi_double_precision,state_vector&
          &,state_dim,mpi_double_precision,cpl_root,cpl_mpi_comms(j),mpi_err)
     !get the tag from the da code
     call mpi_bcast(tag,1,mpi_integer,cpl_root,cpl_mpi_comms(j),mpi_err)
     print*,'Received tag = ',tag
  end do
  print*,'state = ',state_vector


  do i = 1,total_timesteps
     print*,'Timestep = ',i

     do j = 1,ensemble_members_per_proc
        !send the state to da code with mpi_gatherv
        call mpi_gatherv(state_vector,state_dim,mpi_double_precision&
             &,state_vector&
             &,state_dim,state_dim,mpi_double_precision,cpl_root&
             &,cpl_mpi_comms(j),mpi_err)
     end do

     do j = 1,ensemble_members_per_proc
        !get the state back from da code with mpi_gatherv
        call mpi_scatterv(send_null,0,0,mpi_double_precision,state_vector&
             &,state_dim,mpi_double_precision,cpl_root,cpl_mpi_comms(j),mpi_err)
        !get the tag from the da code
        call mpi_bcast(tag,1,mpi_integer,cpl_root,cpl_mpi_comms(j),mpi_err)
        print*,'Received tag = ',tag
        print*,'state = ',state_vector
     end do
  end do
  


















  call mpi_finalize(mpi_err)
  print*,'MINIMAL_MODEL_v5 finished nicely'
  print*,'my communicators were', cpl_mpi_comms(:)
end program minimal_model_v5