[cig-commits] r20391 - seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D

danielpeter at geodynamics.org danielpeter at geodynamics.org
Tue Jun 19 15:21:53 PDT 2012


Author: danielpeter
Date: 2012-06-19 15:21:52 -0700 (Tue, 19 Jun 2012)
New Revision: 20391

Added:
   seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_MPI_interfaces.f90
   seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_inner_outer.f90
Log:
adds setup_inner_outer.f90 and setup_MPI_interfaces.f90

Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_MPI_interfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_MPI_interfaces.f90	                        (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_MPI_interfaces.f90	2012-06-19 22:21:52 UTC (rev 20391)
@@ -0,0 +1,464 @@
+!=====================================================================
+!
+!          S p e c f e m 3 D  G l o b e  V e r s i o n  5 . 1
+!          --------------------------------------------------
+!
+!          Main authors: Dimitri Komatitsch and Jeroen Tromp
+!                        Princeton University, USA
+!             and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+!                            April 2011
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+
+  subroutine setup_MPI_interfaces()
+
+  use meshfem3D_par
+  use create_MPI_interfaces_par
+  implicit none
+
+!  include 'mpif.h'
+
+  ! local parameters
+  integer :: ier,ndim_assemble
+
+  ! temporary buffers for send and receive between faces of the slices and the chunks
+  real(kind=CUSTOM_REAL), dimension(npoin2D_max_all_CM_IC) ::  &
+    buffer_send_faces_scalar,buffer_received_faces_scalar
+
+  ! assigns initial maximum arrays
+  ! for global slices, maximum number of neighbor is around 17 ( 8 horizontal, max of 8 on bottom )
+  integer :: MAX_NEIGHBOURS
+  integer, dimension(:),allocatable :: my_neighbours,nibool_neighbours
+  integer, dimension(:,:),allocatable :: ibool_neighbours
+  integer :: max_nibool
+  real(kind=CUSTOM_REAL),dimension(:),allocatable :: test_flag
+  integer,dimension(:),allocatable :: dummy_i
+  integer :: i,j,k,ispec,iglob
+  ! debug
+  character(len=150) :: filename
+  logical,parameter :: DEBUG_INTERFACES = .false.
+
+  ! estimates a maximum size of needed arrays
+  MAX_NEIGHBOURS = 8 + NCORNERSCHUNKS
+  if( INCLUDE_CENTRAL_CUBE ) MAX_NEIGHBOURS = MAX_NEIGHBOURS + NUMMSGS_FACES
+
+  allocate(my_neighbours(MAX_NEIGHBOURS), &
+          nibool_neighbours(MAX_NEIGHBOURS),stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating my_neighbours array')
+
+  ! estimates initial maximum ibool array
+  max_nibool = npoin2D_max_all_CM_IC * NUMFACES_SHARED &
+               + non_zero_nb_msgs_theor_in_cube*npoin2D_cube_from_slices
+
+  allocate(ibool_neighbours(max_nibool,MAX_NEIGHBOURS), stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating ibool_neighbours')
+
+
+! sets up MPI interfaces
+! crust mantle region
+  if( myrank == 0 ) write(IMAIN,*) 'crust mantle mpi:'
+  allocate(test_flag(NGLOB_CRUST_MANTLE), &
+          stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating test_flag')
+
+  ! sets flag to rank id (+1 to avoid problems with zero rank)
+  test_flag(:) = myrank + 1.0
+
+  ! assembles values
+  call assemble_MPI_scalar_block(myrank,test_flag, &
+            NGLOB_CRUST_MANTLE, &
+            iproc_xi,iproc_eta,ichunk,addressing, &
+            iboolleft_xi_crust_mantle,iboolright_xi_crust_mantle,iboolleft_eta_crust_mantle,iboolright_eta_crust_mantle, &
+            npoin2D_faces_crust_mantle,npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle, &
+            iboolfaces_crust_mantle,iboolcorner_crust_mantle, &
+            iprocfrom_faces,iprocto_faces,imsg_type, &
+            iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners, &
+            buffer_send_faces_scalar,buffer_received_faces_scalar,npoin2D_max_all_CM_IC, &
+            buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar, &
+            NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS, &
+            NPROC_XI,NPROC_ETA,NGLOB1D_RADIAL(IREGION_CRUST_MANTLE), &
+            NGLOB2DMAX_XMIN_XMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_YMIN_YMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_XY,NCHUNKS)
+
+  ! removes own myrank id (+1)
+  test_flag(:) = test_flag(:) - ( myrank + 1.0)
+
+  ! debug: saves array
+  !write(filename,'(a,i6.6)') trim(OUTPUT_FILES)//'/MPI_test_flag_crust_mantle_proc',myrank
+  !call write_VTK_glob_points(NGLOB_CRUST_MANTLE, &
+  !                      xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
+  !                      test_flag,filename)
+
+  allocate(dummy_i(NSPEC_CRUST_MANTLE),stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating dummy_i')
+
+  ! determines neighbor rank for shared faces
+  call get_MPI_interfaces(myrank,NGLOB_CRUST_MANTLE,NSPEC_CRUST_MANTLE, &
+                            test_flag,my_neighbours,nibool_neighbours,ibool_neighbours, &
+                            num_interfaces_crust_mantle,max_nibool_interfaces_crust_mantle, &
+                            max_nibool,MAX_NEIGHBOURS, &
+                            ibool_crust_mantle,&
+                            is_on_a_slice_edge_crust_mantle, &
+                            IREGION_CRUST_MANTLE,.false.,dummy_i,INCLUDE_CENTRAL_CUBE, &
+                            xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle,NPROCTOT)
+
+  deallocate(test_flag)
+  deallocate(dummy_i)
+
+  ! stores MPI interfaces informations
+  allocate(my_neighbours_crust_mantle(num_interfaces_crust_mantle), &
+          nibool_interfaces_crust_mantle(num_interfaces_crust_mantle), &
+          stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating array my_neighbours_crust_mantle etc.')
+  my_neighbours_crust_mantle = -1
+  nibool_interfaces_crust_mantle = 0
+
+  ! copies interfaces arrays
+  if( num_interfaces_crust_mantle > 0 ) then
+    allocate(ibool_interfaces_crust_mantle(max_nibool_interfaces_crust_mantle,num_interfaces_crust_mantle), &
+           stat=ier)
+    if( ier /= 0 ) call exit_mpi(myrank,'error allocating array ibool_interfaces_crust_mantle')
+    ibool_interfaces_crust_mantle = 0
+
+    ! ranks of neighbour processes
+    my_neighbours_crust_mantle(:) = my_neighbours(1:num_interfaces_crust_mantle)
+    ! number of global ibool entries on each interface
+    nibool_interfaces_crust_mantle(:) = nibool_neighbours(1:num_interfaces_crust_mantle)
+    ! global iglob point ids on each interface
+    ibool_interfaces_crust_mantle(:,:) = ibool_neighbours(1:max_nibool_interfaces_crust_mantle,1:num_interfaces_crust_mantle)
+  else
+    ! dummy allocation (fortran90 should allow allocate statement with zero array size)
+    max_nibool_interfaces_crust_mantle = 0
+    allocate(ibool_interfaces_crust_mantle(0,0),stat=ier)
+  endif
+
+  ! debug: outputs MPI interface
+  if( DEBUG_INTERFACES ) then
+  do i=1,num_interfaces_crust_mantle
+    write(filename,'(a,i6.6,a,i2.2)') trim(OUTPUT_FILES)//'/MPI_points_crust_mantle_proc',myrank, &
+                    '_',my_neighbours_crust_mantle(i)
+    call write_VTK_data_points(NGLOB_crust_mantle, &
+                      xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
+                      ibool_interfaces_crust_mantle(1:nibool_interfaces_crust_mantle(i),i), &
+                      nibool_interfaces_crust_mantle(i),filename)
+  enddo
+  call sync_all()
+  endif
+
+  ! checks addressing
+  call test_MPI_neighbours(IREGION_CRUST_MANTLE, &
+                              num_interfaces_crust_mantle,max_nibool_interfaces_crust_mantle, &
+                              my_neighbours_crust_mantle,nibool_interfaces_crust_mantle, &
+                              ibool_interfaces_crust_mantle)
+
+  ! allocates MPI buffers
+  ! crust mantle
+  allocate(buffer_send_vector_crust_mantle(NDIM,max_nibool_interfaces_crust_mantle,num_interfaces_crust_mantle), &
+          buffer_recv_vector_crust_mantle(NDIM,max_nibool_interfaces_crust_mantle,num_interfaces_crust_mantle), &
+          request_send_vector_crust_mantle(num_interfaces_crust_mantle), &
+          request_recv_vector_crust_mantle(num_interfaces_crust_mantle), &
+          stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating array buffer_send_vector_crust_mantle etc.')
+
+  ! checks with assembly of test fields
+  call test_MPI_cm()
+
+
+! outer core region
+  if( myrank == 0 ) write(IMAIN,*) 'outer core mpi:'
+
+  allocate(test_flag(NGLOB_OUTER_CORE), &
+          stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating test_flag outer core')
+
+  ! sets flag to rank id (+1 to avoid problems with zero rank)
+  test_flag(:) = myrank + 1.0
+
+  ! assembles values
+  call assemble_MPI_scalar_block(myrank,test_flag, &
+            NGLOB_OUTER_CORE, &
+            iproc_xi,iproc_eta,ichunk,addressing, &
+            iboolleft_xi_outer_core,iboolright_xi_outer_core,iboolleft_eta_outer_core,iboolright_eta_outer_core, &
+            npoin2D_faces_outer_core,npoin2D_xi_outer_core,npoin2D_eta_outer_core, &
+            iboolfaces_outer_core,iboolcorner_outer_core, &
+            iprocfrom_faces,iprocto_faces,imsg_type, &
+            iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners, &
+            buffer_send_faces_scalar,buffer_received_faces_scalar,npoin2D_max_all_CM_IC, &
+            buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar, &
+            NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS, &
+            NPROC_XI,NPROC_ETA,NGLOB1D_RADIAL(IREGION_OUTER_CORE), &
+            NGLOB2DMAX_XMIN_XMAX(IREGION_OUTER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_OUTER_CORE),NGLOB2DMAX_XY,NCHUNKS)
+
+
+  ! removes own myrank id (+1)
+  test_flag(:) = test_flag(:) - ( myrank + 1.0)
+
+  ! debug: saves array
+  !write(filename,'(a,i6.6)') trim(OUTPUT_FILES)//'/MPI_test_flag_outer_core_proc',myrank
+  !call write_VTK_glob_points(NGLOB_OUTER_CORE, &
+  !                      xstore_outer_core,ystore_outer_core,zstore_outer_core, &
+  !                      test_flag,filename)
+
+  allocate(dummy_i(NSPEC_OUTER_CORE),stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating dummy_i')
+
+  ! determines neighbor rank for shared faces
+  call get_MPI_interfaces(myrank,NGLOB_OUTER_CORE,NSPEC_OUTER_CORE, &
+                            test_flag,my_neighbours,nibool_neighbours,ibool_neighbours, &
+                            num_interfaces_outer_core,max_nibool_interfaces_outer_core, &
+                            max_nibool,MAX_NEIGHBOURS, &
+                            ibool_outer_core,&
+                            is_on_a_slice_edge_outer_core, &
+                            IREGION_OUTER_CORE,.false.,dummy_i,INCLUDE_CENTRAL_CUBE, &
+                            xstore_outer_core,ystore_outer_core,zstore_outer_core,NPROCTOT)
+
+  deallocate(test_flag)
+  deallocate(dummy_i)
+
+  ! stores MPI interfaces informations
+  allocate(my_neighbours_outer_core(num_interfaces_outer_core), &
+          nibool_interfaces_outer_core(num_interfaces_outer_core), &
+          stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating array my_neighbours_outer_core etc.')
+  my_neighbours_outer_core = -1
+  nibool_interfaces_outer_core = 0
+
+  ! copies interfaces arrays
+  if( num_interfaces_outer_core > 0 ) then
+    allocate(ibool_interfaces_outer_core(max_nibool_interfaces_outer_core,num_interfaces_outer_core), &
+           stat=ier)
+    if( ier /= 0 ) call exit_mpi(myrank,'error allocating array ibool_interfaces_outer_core')
+    ibool_interfaces_outer_core = 0
+
+    ! ranks of neighbour processes
+    my_neighbours_outer_core(:) = my_neighbours(1:num_interfaces_outer_core)
+    ! number of global ibool entries on each interface
+    nibool_interfaces_outer_core(:) = nibool_neighbours(1:num_interfaces_outer_core)
+    ! global iglob point ids on each interface
+    ibool_interfaces_outer_core(:,:) = ibool_neighbours(1:max_nibool_interfaces_outer_core,1:num_interfaces_outer_core)
+  else
+    ! dummy allocation (fortran90 should allow allocate statement with zero array size)
+    max_nibool_interfaces_outer_core = 0
+    allocate(ibool_interfaces_outer_core(0,0),stat=ier)
+  endif
+
+  ! debug: outputs MPI interface
+  if( DEBUG_INTERFACES ) then
+  do i=1,num_interfaces_outer_core
+    write(filename,'(a,i6.6,a,i2.2)') trim(OUTPUT_FILES)//'/MPI_points_outer_core_proc',myrank, &
+                    '_',my_neighbours_outer_core(i)
+    call write_VTK_data_points(NGLOB_OUTER_CORE, &
+                      xstore_outer_core,ystore_outer_core,zstore_outer_core, &
+                      ibool_interfaces_outer_core(1:nibool_interfaces_outer_core(i),i), &
+                      nibool_interfaces_outer_core(i),filename)
+  enddo
+  call sync_all()
+  endif
+
+  ! checks addressing
+  call test_MPI_neighbours(IREGION_OUTER_CORE, &
+                              num_interfaces_outer_core,max_nibool_interfaces_outer_core, &
+                              my_neighbours_outer_core,nibool_interfaces_outer_core, &
+                              ibool_interfaces_outer_core)
+
+  ! allocates MPI buffers
+  ! outer core
+  allocate(buffer_send_scalar_outer_core(max_nibool_interfaces_outer_core,num_interfaces_outer_core), &
+          buffer_recv_scalar_outer_core(max_nibool_interfaces_outer_core,num_interfaces_outer_core), &
+          request_send_scalar_outer_core(num_interfaces_outer_core), &
+          request_recv_scalar_outer_core(num_interfaces_outer_core), &
+          stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating array buffer_send_vector_outer_core etc.')
+
+  ! checks with assembly of test fields
+  call test_MPI_oc()
+
+
+! inner core
+  if( myrank == 0 ) write(IMAIN,*) 'inner core mpi:'
+
+  allocate(test_flag(NGLOB_INNER_CORE), &
+          stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating test_flag inner core')
+
+  ! sets flag to rank id (+1 to avoid problems with zero rank)
+  test_flag(:) = 0.0
+  do ispec=1,NSPEC_INNER_CORE
+    ! suppress fictitious elements in central cube
+    if(idoubling_inner_core(ispec) == IFLAG_IN_FICTITIOUS_CUBE) cycle
+    ! sets flags
+    do k = 1,NGLLZ
+      do j = 1,NGLLY
+        do i = 1,NGLLX
+          iglob = ibool_inner_core(i,j,k,ispec)
+          test_flag(iglob) = myrank + 1.0
+        enddo
+      enddo
+    enddo
+  enddo
+
+  ! assembles values
+  call assemble_MPI_scalar_block(myrank,test_flag, &
+            NGLOB_INNER_CORE, &
+            iproc_xi,iproc_eta,ichunk,addressing, &
+            iboolleft_xi_inner_core,iboolright_xi_inner_core,iboolleft_eta_inner_core,iboolright_eta_inner_core, &
+            npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core, &
+            iboolfaces_inner_core,iboolcorner_inner_core, &
+            iprocfrom_faces,iprocto_faces,imsg_type, &
+            iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners, &
+            buffer_send_faces_scalar,buffer_received_faces_scalar,npoin2D_max_all_CM_IC, &
+            buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar, &
+            NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS, &
+            NPROC_XI,NPROC_ETA,NGLOB1D_RADIAL(IREGION_INNER_CORE), &
+            NGLOB2DMAX_XMIN_XMAX(IREGION_INNER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_INNER_CORE),NGLOB2DMAX_XY,NCHUNKS)
+
+  ! debug: saves array
+  !write(filename,'(a,i6.6)') trim(OUTPUT_FILES)//'/MPI_test_flag_inner_core_A_proc',myrank
+  !call write_VTK_glob_points(NGLOB_INNER_CORE, &
+  !                      xstore_inner_core,ystore_inner_core,zstore_inner_core, &
+  !                      test_flag,filename)
+
+  ! debug: idoubling inner core
+  if( DEBUG_INTERFACES ) then
+    write(filename,'(a,i6.6)') trim(OUTPUT_FILES)//'/MPI_idoubling_inner_core_proc',myrank
+    call write_VTK_data_elem_i(NSPEC_INNER_CORE,NGLOB_INNER_CORE, &
+                            xstore_inner_core,ystore_inner_core,zstore_inner_core, &
+                            ibool_inner_core, &
+                            idoubling_inner_core,filename)
+    call sync_all()
+  endif
+
+  ! including central cube
+  if(INCLUDE_CENTRAL_CUBE) then
+    ! user output
+    if( myrank == 0 ) write(IMAIN,*) 'inner core with central cube mpi:'
+
+    ! test_flag is a scalar, not a vector
+    ndim_assemble = 1
+
+    ! use central cube buffers to assemble the inner core mass matrix with the central cube
+    call assemble_MPI_central_cube_block(ichunk,nb_msgs_theor_in_cube, sender_from_slices_to_cube, &
+                 npoin2D_cube_from_slices, buffer_all_cube_from_slices, &
+                 buffer_slices, buffer_slices2, ibool_central_cube, &
+                 receiver_cube_from_slices, ibool_inner_core, &
+                 idoubling_inner_core, NSPEC_INNER_CORE, &
+                 ibelm_bottom_inner_core, NSPEC2D_BOTTOM(IREGION_INNER_CORE), &
+                 NGLOB_INNER_CORE, &
+                 test_flag,ndim_assemble, &
+                 iproc_eta,addressing,NCHUNKS,NPROC_XI,NPROC_ETA)
+  endif
+
+  ! removes own myrank id (+1)
+  test_flag = test_flag - ( myrank + 1.0)
+  where( test_flag < 0.0 ) test_flag = 0.0
+
+  ! debug: saves array
+  !write(filename,'(a,i6.6)') trim(OUTPUT_FILES)//'/MPI_test_flag_inner_core_B_proc',myrank
+  !call write_VTK_glob_points(NGLOB_INNER_CORE, &
+  !                    xstore_inner_core,ystore_inner_core,zstore_inner_core, &
+  !                    test_flag,filename)
+  !call sync_all()
+
+  ! in sequential order, for testing purpose
+  do i=0,NPROCTOT - 1
+    if( myrank == i ) then
+      ! gets new interfaces for inner_core without central cube yet
+      ! determines neighbor rank for shared faces
+      call get_MPI_interfaces(myrank,NGLOB_INNER_CORE,NSPEC_INNER_CORE, &
+                            test_flag,my_neighbours,nibool_neighbours,ibool_neighbours, &
+                            num_interfaces_inner_core,max_nibool_interfaces_inner_core, &
+                            max_nibool,MAX_NEIGHBOURS, &
+                            ibool_inner_core,&
+                            is_on_a_slice_edge_inner_core, &
+                            IREGION_INNER_CORE,.false.,idoubling_inner_core,INCLUDE_CENTRAL_CUBE, &
+                            xstore_inner_core,ystore_inner_core,zstore_inner_core,NPROCTOT)
+
+    endif
+    call sync_all()
+  enddo
+
+
+  deallocate(test_flag)
+  call sync_all()
+
+  ! stores MPI interfaces informations
+  allocate(my_neighbours_inner_core(num_interfaces_inner_core), &
+          nibool_interfaces_inner_core(num_interfaces_inner_core), &
+          stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating array my_neighbours_inner_core etc.')
+  my_neighbours_inner_core = -1
+  nibool_interfaces_inner_core = 0
+
+  ! copies interfaces arrays
+  if( num_interfaces_inner_core > 0 ) then
+    allocate(ibool_interfaces_inner_core(max_nibool_interfaces_inner_core,num_interfaces_inner_core), &
+           stat=ier)
+    if( ier /= 0 ) call exit_mpi(myrank,'error allocating array ibool_interfaces_inner_core')
+    ibool_interfaces_inner_core = 0
+
+    ! ranks of neighbour processes
+    my_neighbours_inner_core(:) = my_neighbours(1:num_interfaces_inner_core)
+    ! number of global ibool entries on each interface
+    nibool_interfaces_inner_core(:) = nibool_neighbours(1:num_interfaces_inner_core)
+    ! global iglob point ids on each interface
+    ibool_interfaces_inner_core(:,:) = ibool_neighbours(1:max_nibool_interfaces_inner_core,1:num_interfaces_inner_core)
+  else
+    ! dummy allocation (fortran90 should allow allocate statement with zero array size)
+    max_nibool_interfaces_inner_core = 0
+    allocate(ibool_interfaces_inner_core(0,0),stat=ier)
+  endif
+
+  ! debug: saves MPI interfaces
+  if( DEBUG_INTERFACES ) then
+  do i=1,num_interfaces_inner_core
+    write(filename,'(a,i6.6,a,i2.2)') trim(OUTPUT_FILES)//'/MPI_points_inner_core_proc',myrank, &
+                    '_',my_neighbours_inner_core(i)
+    call write_VTK_data_points(NGLOB_INNER_CORE, &
+                      xstore_inner_core,ystore_inner_core,zstore_inner_core, &
+                      ibool_interfaces_inner_core(1:nibool_interfaces_inner_core(i),i), &
+                      nibool_interfaces_inner_core(i),filename)
+  enddo
+  call sync_all()
+  endif
+
+  ! checks addressing
+  call test_MPI_neighbours(IREGION_INNER_CORE, &
+                              num_interfaces_inner_core,max_nibool_interfaces_inner_core, &
+                              my_neighbours_inner_core,nibool_interfaces_inner_core, &
+                              ibool_interfaces_inner_core)
+
+  ! allocates MPI buffers
+  ! inner core
+  allocate(buffer_send_vector_inner_core(NDIM,max_nibool_interfaces_inner_core,num_interfaces_inner_core), &
+          buffer_recv_vector_inner_core(NDIM,max_nibool_interfaces_inner_core,num_interfaces_inner_core), &
+          request_send_vector_inner_core(num_interfaces_inner_core), &
+          request_recv_vector_inner_core(num_interfaces_inner_core), &
+          stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating array buffer_send_vector_inner_core etc.')
+
+  ! checks with assembly of test fields
+  call test_MPI_ic()
+
+  ! synchronizes MPI processes
+  call sync_all()
+
+  ! frees temporary array
+  deallocate(ibool_neighbours)
+  deallocate(my_neighbours,nibool_neighbours)
+
+  end subroutine setup_MPI_interfaces

Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_inner_outer.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_inner_outer.f90	                        (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_inner_outer.f90	2012-06-19 22:21:52 UTC (rev 20391)
@@ -0,0 +1,165 @@
+!=====================================================================
+!
+!          S p e c f e m 3 D  G l o b e  V e r s i o n  5 . 1
+!          --------------------------------------------------
+!
+!          Main authors: Dimitri Komatitsch and Jeroen Tromp
+!                        Princeton University, USA
+!             and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+!                            April 2011
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+
+  subroutine setup_Inner_Outer()
+
+  use meshfem3D_par
+  use create_MPI_interfaces_par
+  implicit none
+
+  ! local parameters
+  real :: percentage_edge
+  integer :: ier,ispec,iinner,iouter
+  ! debug
+  character(len=150) :: filename
+  logical,parameter :: DEBUG_INTERFACES = .false.
+
+  ! stores inner / outer elements
+  !
+  ! note: arrays is_on_a_slice_edge_.. have flags set for elements which need to
+  !         communicate with other MPI processes
+
+  ! crust_mantle
+  nspec_outer_crust_mantle = count( is_on_a_slice_edge_crust_mantle )
+  nspec_inner_crust_mantle = NSPEC_CRUST_MANTLE - nspec_outer_crust_mantle
+
+  num_phase_ispec_crust_mantle = max(nspec_inner_crust_mantle,nspec_outer_crust_mantle)
+
+  allocate(phase_ispec_inner_crust_mantle(num_phase_ispec_crust_mantle,2),stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating array phase_ispec_inner_crust_mantle')
+
+  phase_ispec_inner_crust_mantle(:,:) = 0
+  iinner = 0
+  iouter = 0
+  do ispec=1,NSPEC_CRUST_MANTLE
+    if( is_on_a_slice_edge_crust_mantle(ispec) ) then
+      ! outer element
+      iouter = iouter + 1
+      phase_ispec_inner_crust_mantle(iouter,1) = ispec
+    else
+      ! inner element
+      iinner = iinner + 1
+      phase_ispec_inner_crust_mantle(iinner,2) = ispec
+    endif
+  enddo
+
+  ! outer_core
+  nspec_outer_outer_core = count( is_on_a_slice_edge_outer_core )
+  nspec_inner_outer_core = NSPEC_OUTER_CORE - nspec_outer_outer_core
+
+  num_phase_ispec_outer_core = max(nspec_inner_outer_core,nspec_outer_outer_core)
+
+  allocate(phase_ispec_inner_outer_core(num_phase_ispec_outer_core,2),stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating array phase_ispec_inner_outer_core')
+
+  phase_ispec_inner_outer_core(:,:) = 0
+  iinner = 0
+  iouter = 0
+  do ispec=1,NSPEC_OUTER_CORE
+    if( is_on_a_slice_edge_outer_core(ispec) ) then
+      ! outer element
+      iouter = iouter + 1
+      phase_ispec_inner_outer_core(iouter,1) = ispec
+    else
+      ! inner element
+      iinner = iinner + 1
+      phase_ispec_inner_outer_core(iinner,2) = ispec
+    endif
+  enddo
+
+  ! inner_core
+  nspec_outer_inner_core = count( is_on_a_slice_edge_inner_core )
+  nspec_inner_inner_core = NSPEC_INNER_CORE - nspec_outer_inner_core
+
+  num_phase_ispec_inner_core = max(nspec_inner_inner_core,nspec_outer_inner_core)
+
+  allocate(phase_ispec_inner_inner_core(num_phase_ispec_inner_core,2),stat=ier)
+  if( ier /= 0 ) call exit_mpi(myrank,'error allocating array phase_ispec_inner_inner_core')
+
+  phase_ispec_inner_inner_core(:,:) = 0
+  iinner = 0
+  iouter = 0
+  do ispec=1,NSPEC_INNER_CORE
+    if( is_on_a_slice_edge_inner_core(ispec) ) then
+      ! outer element
+      iouter = iouter + 1
+      phase_ispec_inner_inner_core(iouter,1) = ispec
+    else
+      ! inner element
+      iinner = iinner + 1
+      phase_ispec_inner_inner_core(iinner,2) = ispec
+    endif
+  enddo
+
+  ! user output
+  if(myrank == 0) then
+
+    write(IMAIN,*)
+    write(IMAIN,*) 'for overlapping of communications with calculations:'
+    write(IMAIN,*)
+
+    percentage_edge = 100. * nspec_outer_crust_mantle / real(NSPEC_CRUST_MANTLE)
+    write(IMAIN,*) 'percentage of edge elements in crust/mantle ',percentage_edge,'%'
+    write(IMAIN,*) 'percentage of volume elements in crust/mantle ',100. - percentage_edge,'%'
+    write(IMAIN,*)
+
+    percentage_edge = 100.* nspec_outer_outer_core / real(NSPEC_OUTER_CORE)
+    write(IMAIN,*) 'percentage of edge elements in outer core ',percentage_edge,'%'
+    write(IMAIN,*) 'percentage of volume elements in outer core ',100. - percentage_edge,'%'
+    write(IMAIN,*)
+
+    percentage_edge = 100. * nspec_outer_inner_core / real(NSPEC_INNER_CORE)
+    write(IMAIN,*) 'percentage of edge elements in inner core ',percentage_edge,'%'
+    write(IMAIN,*) 'percentage of volume elements in inner core ',100. - percentage_edge,'%'
+    write(IMAIN,*)
+
+  endif
+
+  ! debug: saves element flags
+  if( DEBUG_INTERFACES ) then
+    ! crust mantle
+    write(filename,'(a,i6.6)') trim(OUTPUT_FILES)//'/MPI_innerouter_crust_mantle_proc',myrank
+    call write_VTK_data_elem_l(NSPEC_CRUST_MANTLE,NGLOB_CRUST_MANTLE, &
+                              xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
+                              ibool_crust_mantle, &
+                              is_on_a_slice_edge_crust_mantle,filename)
+    ! outer core
+    write(filename,'(a,i6.6)') trim(OUTPUT_FILES)//'/MPI_innerouter_outer_core_proc',myrank
+    call write_VTK_data_elem_l(NSPEC_OUTER_CORE,NGLOB_OUTER_CORE, &
+                              xstore_outer_core,ystore_outer_core,zstore_outer_core, &
+                              ibool_outer_core, &
+                              is_on_a_slice_edge_outer_core,filename)
+    ! inner core
+    write(filename,'(a,i6.6)') trim(OUTPUT_FILES)//'/MPI_innerouter_inner_core_proc',myrank
+    call write_VTK_data_elem_l(NSPEC_INNER_CORE,NGLOB_INNER_CORE, &
+                              xstore_inner_core,ystore_inner_core,zstore_inner_core, &
+                              ibool_inner_core, &
+                              is_on_a_slice_edge_inner_core,filename)
+  endif
+
+  end subroutine setup_Inner_Outer



More information about the CIG-COMMITS mailing list