[cig-commits] r20216 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src

[surendra at geodynamics.org]

Author: surendra
Date: 2012-05-28 22:59:26 -0700 (Mon, 28 May 2012)
New Revision: 20216

Modified:
   seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Parallelized fault time-series output

Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90	2012-05-29 02:10:05 UTC (rev 20215)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90	2012-05-29 05:59:26 UTC (rev 20216)
@@ -39,8 +39,12 @@
   ! outputs on selected fault nodes at every time step:
   ! slip, slip velocity, fault stresses
   type dataT_type
-    integer                                    :: npoin
+    integer                                    :: npoin,npoin_local
     integer, dimension(:), pointer             :: iglob   ! on-fault global index of output nodes
+    integer, dimension(:,:), pointer             :: iglob_all
+    integer, dimension(:), pointer             :: islice,glob_indx
+    real(kind=CUSTOM_REAL), dimension(:), pointer  :: dist
+    real(kind=CUSTOM_REAL), dimension(:,:), pointer  :: dist_all
     real(kind=CUSTOM_REAL), dimension(:,:), pointer  :: d1,v1,t1,d2,v2,t2,t3,theta
     character(len=70), dimension(:), pointer   :: name
   end type dataT_type
@@ -958,7 +962,7 @@
   ! write dataXZ every NSNAP time step
   if ( mod(it,NSNAP) == 0) then
     if(.NOT. PARALLEL_FAULT) then
-      call write_dataXZ(bc%dataXZ,it,iflt)
+      if(bc%nspec > 0) call write_dataXZ(bc%dataXZ,it,iflt)
     else
       call gather_dataXZ(bc)
       if(myrank==0) call write_dataXZ(bc%dataXZ_all,it,iflt)
@@ -1166,6 +1170,7 @@
 !===============================================================
 ! OUTPUTS
 subroutine init_dataT(DataT,coord,nglob,NT,iflt)
+  use specfem_par, only : NPROC,myrank
   ! NT = total number of time steps
 
   integer, intent(in) :: nglob,NT,iflt
@@ -1176,6 +1181,8 @@
   integer :: i, iglob , IIN, ier, jflt, np, k
   character(len=70) :: tmpname
 
+  integer :: ipoin, ipoin_local, iproc
+
   !  1. read fault output coordinates from user file, 
   !  2. define iglob: the fault global index of the node nearest to user
   !     requested coordinate
@@ -1194,6 +1201,7 @@
 
   allocate(DataT%iglob(DataT%npoin))
   allocate(DataT%name(DataT%npoin))
+  allocate(DataT%dist(DataT%npoin)) !Surendra : for parallel fault
 
   open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
   if( ier /= 0 ) stop 'error opening FAULT_STATIONS file'
@@ -1216,8 +1224,43 @@
         DataT%iglob(k) = iglob   
       endif
     enddo
+    DataT%dist(k) = distkeep !Surendra : for parallel fault
   enddo
 
+ !Surendra : for Parallel fault
+  if(PARALLEL_FAULT) then
+    allocate(DataT%islice(DataT%npoin)) 
+    allocate(DataT%iglob_all(DataT%npoin,0:NPROC-1))
+    allocate(DataT%dist_all(DataT%npoin,0:NPROC-1))
+    call gather_all_i(DataT%iglob,DataT%npoin,DataT%iglob_all,DataT%npoin,NPROC)
+    call gather_all_cr(DataT%dist,DataT%npoin,DataT%dist_all,DataT%npoin,NPROC)
+    if(myrank==0) then
+      do ipoin = 1,DataT%npoin
+        iproc = minloc(DataT%dist_all(ipoin,:), 1) - 1
+        DataT%islice(ipoin) = iproc
+        DataT%iglob(ipoin) = DataT%iglob_all(ipoin,iproc)
+      enddo
+    endif
+
+    call bcast_all_i(DataT%islice,DataT%npoin)
+    call bcast_all_i(DataT%iglob,DataT%npoin)
+
+    DataT%npoin_local = 0
+    do ipoin = 1,DataT%npoin
+      if(myrank == DataT%islice(ipoin)) DataT%npoin_local = DataT%npoin_local + 1
+    enddo
+    allocate(DataT%glob_indx(DataT%npoin_local)) 
+    do ipoin = 1,DataT%npoin
+      if(myrank == DataT%islice(ipoin)) then
+        ipoin_local = ipoin_local + 1
+        DataT%glob_indx(ipoin_local) = ipoin
+      endif
+    enddo
+  else
+    DataT%npoin_local = DataT%npoin
+  endif !Parallel_fault
+
+
   !  3. allocate arrays and set to zero
   allocate(DataT%d1(NT,DataT%npoin))
   allocate(DataT%v1(NT,DataT%npoin))
@@ -1248,9 +1291,14 @@
   real(kind=CUSTOM_REAL), dimension(:), intent(in) :: theta
   integer, intent(in) :: itime
 
-  integer :: i,k
+  integer :: i,k,ipoin
 
-  do i=1,dataT%npoin
+  do ipoin=1,dataT%npoin_local
+    if (PARALLEL_FAULT) then
+      i = DataT%glob_indx(ipoin)
+    else
+      i = ipoin
+    endif
     k = dataT%iglob(i)
     dataT%d1(itime,i) = d(1,k)
     dataT%d2(itime,i) = d(2,k)
@@ -1285,7 +1333,7 @@
   real(kind=CUSTOM_REAL), intent(in) :: DT
   integer, intent(in) :: NT
 
-  integer   :: i,k,IOUT
+  integer   :: i,k,IOUT,ipoin
   character :: NTchar*5
   integer ::  today(3), now(3)
 
@@ -1298,7 +1346,12 @@
   NTchar = adjustl(NTchar)
 
 1 format(I5)  
-  do i=1,dataT%npoin
+  do ipoin=1,dataT%npoin_local
+    if (PARALLEL_FAULT) then
+      i = DataT%glob_indx(ipoin)
+    else
+      i = ipoin
+    endif
 
     open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
     write(IOUT,*) "# problem=TPV102"