[cig-commits] r20209 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: decompose_mesh_SCOTCH src
surendra at geodynamics.org
surendra at geodynamics.org
Sun May 27 15:20:54 PDT 2012
Author: surendra
Date: 2012-05-27 15:20:53 -0700 (Sun, 27 May 2012)
New Revision: 20209
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/iterate_time.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/prepare_timerun.f90
Log:
Merged fault parallelization into new_fault_db branch with flag for parallel fault. Default is serial fault
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-05-27 07:40:22 UTC (rev 20208)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-05-27 22:20:53 UTC (rev 20209)
@@ -6,7 +6,7 @@
use fault_scotch
implicit none
-
+
include './scotchf.h'
! number of partitions
@@ -83,6 +83,8 @@
integer :: q_flag,aniso_flag,idomain_id
double precision :: vp,vs,rho
+ logical, parameter :: PARALLEL_FAULT = .false.
+
contains
!----------------------------------------------------------------------------------------------
@@ -389,10 +391,13 @@
if( nspec2D_moho > 0 ) print*, ' nspec2D_moho = ', nspec2D_moho
call read_fault_files(localpath_name)
- if (faults_exist) then
+ if (allocated(faults)) then
call save_nodes_coords(nodes_coords,nnodes)
call close_faults(nodes_coords,elmnts,nspec,nnodes,esize)
-!JPA parallel fault call reorder_fault_elements(nodes_coords,nnodes)
+ !SNS
+ if(PARALLEL_FAULT) then
+ call reorder_fault_elements(nodes_coords,elmnts,nnodes,esize,nspec)
+ endif
end if
end subroutine read_mesh_files
@@ -530,10 +535,13 @@
! count_def_mat, mat(1,:) , mat_prop, &
! sup_neighbour, nsize, &
! nparts, part, nfaces_coupled, faces_coupled)
-
- ! move all fault elements to the same partition (proc=0)
- call fault_repartition (nspec, nnodes, elmnts, nsize, nparts, part, esize)
+ if(.not. PARALLEL_FAULT) then
+ ! move all fault elements to the same partition (proc=0)
+ call fault_repartition (nspec, nnodes, elmnts, nsize, nparts, part, esize)
+ else
+ call fault_repartition_parallel (nspec, nnodes, elmnts, nsize, nparts, part, esize)
+ endif
! re-partitioning puts moho-surface coupled elements into same partition
call moho_surface_repartitioning (nspec, nnodes, elmnts, &
@@ -650,18 +658,18 @@
close(15)
- if (faults_exist) cycle
-!jpa: TO DO: move all this to a subroutine in fault_scotch.f90
+ if (.not. allocated(faults)) cycle
! write fault database
write(prname, "(i6.6,'_Database_fault')") ipart
open(unit=16,file=outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname,&
- status='unknown', action='write', form='formatted', iostat = ierr)
+ status='replace', action='write', form='formatted', iostat = ierr)
if( ierr /= 0 ) then
print*,'error file open:',outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname
print*
print*,'check if path exists:',outputpath_name(1:len_trim(outputpath_name))
stop
endif
+print*,prname
call write_fault_database(16, ipart, nspec, &
glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, part)
@@ -699,3 +707,4 @@
end module decompose_mesh_SCOTCH
+
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2012-05-27 07:40:22 UTC (rev 20208)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2012-05-27 22:20:53 UTC (rev 20209)
@@ -1,8 +1,8 @@
module fault_scotch
+ implicit none
+ include "../constants.h"
- implicit none
private
-
type fault_type
private
integer :: nspec
@@ -10,32 +10,21 @@
integer, dimension(:,:), pointer :: inodes1, inodes2
end type fault_type
- ! Singleton design pattern: only one instance of this object can be created
type(fault_type), allocatable, save :: faults(:)
-
double precision, dimension(:,:), allocatable, save :: nodes_coords_open
+
integer, parameter :: long = SELECTED_INT_KIND(18)
double precision, parameter :: FAULT_GAP_TOLERANCE = 1.0d0
- public :: read_fault_files, fault_repartition, close_faults, write_fault_database, &
- save_nodes_coords, nodes_coords_open, faults_exist
-!JPA parallel fault , reorder_fault_elements
-
+ public :: read_fault_files, fault_repartition, fault_repartition_parallel, close_faults, write_fault_database, &
+ save_nodes_coords, nodes_coords_open, faults, reorder_fault_elements
CONTAINS
-
!==========================================================================================
- logical function faults_exist()
- faults_exist = allocated(faults)
+ Subroutine read_fault_files(localpath_name)
- end function faults_exist
-
-!==========================================================================================
-
- subroutine read_fault_files(localpath_name)
-
character(len=256),intent(in) :: localpath_name
integer :: nbfaults, iflt, ier
@@ -115,12 +104,13 @@
! ---------------------------------------------------------------------------------------------------
-! Saving nodes_coords with open split nodes to be used in SESAME for ibool_fault_side1 and side2
- subroutine save_nodes_coords(nodes_coords)
+! Saving nodes_coords to be used in SESAME for ibool_fault_side1 and side2
+ subroutine save_nodes_coords(nodes_coords,nnodes)
- double precision, dimension(:,:), intent(in) :: nodes_coords
+ integer, intent(in) :: nnodes
+ double precision, dimension(3,nnodes), intent(in) :: nodes_coords
- allocate(nodes_coords_open(size(nodes_coords,1),size(nodes_coords,2)))
+ allocate(nodes_coords_open(3,nnodes))
nodes_coords_open = nodes_coords
end subroutine save_nodes_coords
@@ -213,53 +203,217 @@
end subroutine close_fault_single
!===================================================================================================
-!JPA parallel fault
-! subroutine reorder_fault_elements(nodes_coords,nnodes)
-!!JPA ...
-! do iflt=1,size(faults)
-! call reorder_fault_elements_single(faults(iflt),&
-! elmnts,nodes_coords,nnodes,esize,nelmnts)
-! enddo
-!
-! end subroutine reorder_fault_elements
-!
-!! ---------------------------------------------------------------------------------------------------
-! subroutine reorder_fault_elements_single(f)
-!
-! ! compute element-face centers for fault side 1
-! do i=1,f%nspec
-! e = f%ispec1(i)
-! do k=1,4
-! iglob = f%inodes1(k,e)
-! xyz(:,k) = nodes_coords(:,iglob)
-! enddo
-! xyz_c(:,k) = ! ... average of xyz
-! enddo
-! ! reorder
-! call lex_order(xyz_c,new_index_list)
-! f%ispec1 = f%ispec1(new_index_list)
-!
-! ! repeat for fault side 2
-! ! ....
-!
-! end subroutine reorder_fault_elements_single
-!
-!! ---------------------------------------------------------------------------------------------------
-! subroutine lex_order()
-!! JPA lexicographic ordering inspired from src/get_global.f90
-! end subroutine lex_order
+ subroutine reorder_fault_elements(nodes_coords,elmnts,nnodes,esize,nelmnts)
+ integer ,intent(in) :: nnodes, esize
+ integer(long), intent(in) :: nelmnts
+ integer, dimension(esize,nelmnts), intent(in) :: elmnts
+ double precision,dimension(3,nnodes), intent(inout) :: nodes_coords
+ integer :: iflt
+
+ do iflt=1,size(faults)
+ call reorder_fault_elements_single(faults(iflt),&
+ elmnts,nodes_coords,nnodes,esize,nelmnts,faults(iflt)%nspec)
+ enddo
+
+ end subroutine reorder_fault_elements
+
+! ---------------------------------------------------------------------------------------------------
+ subroutine reorder_fault_elements_single(f,elmnts,nodes_coords,nnodes,esize,nelmnts,nspec)
+
+ integer ,intent(in) :: nnodes, esize
+ integer(long), intent(in) :: nelmnts
+ integer, dimension(esize,nelmnts), intent(in) :: elmnts
+ double precision,dimension(3,nnodes), intent(inout) :: nodes_coords
+ type(fault_type), intent(inout) :: f
+
+ double precision, dimension(3,4) :: xyz
+ double precision, dimension(:,:), allocatable :: xyz_c
+
+ double precision :: dist
+ integer :: iglob1, iglob2, i, j, k
+ logical :: found_it
+
+ integer :: e,iglob
+
+ integer :: nspec
+
+ integer, dimension(nspec) :: new_index_list
+
+ allocate(xyz_c(3,f%nspec))
+ xyz_c=0._CUSTOM_REAL
+
+ ! compute element-face centers for fault side 1
+ do e=1,f%nspec
+ do k=1,4
+ iglob = f%inodes1(k,e)
+ xyz(:,k) = nodes_coords(:,iglob)
+ enddo
+ xyz_c(1,e) = sum(xyz(1,:))/4
+ xyz_c(2,e) = sum(xyz(2,:))/4
+ xyz_c(3,e) = sum(xyz(3,:))/4
+ enddo
+ ! reorder
+ call lex_order(xyz_c,new_index_list,nnodes,f%nspec)
+ f%ispec1 = f%ispec1(new_index_list)
+ f%inodes1 = f%inodes1(:,new_index_list)
+
+ ! repeat for fault side 2
+ do e=1,f%nspec
+ do k=1,4
+ iglob = f%inodes2(k,e)
+ xyz(:,k) = nodes_coords(:,iglob)
+ enddo
+ xyz_c(1,e) = sum(xyz(1,:))/4
+ xyz_c(2,e) = sum(xyz(2,:))/4
+ xyz_c(3,e) = sum(xyz(3,:))/4
+ enddo
+ ! reorder
+ call lex_order(xyz_c,new_index_list,nnodes,f%nspec)
+ f%ispec2 = f%ispec2(new_index_list)
+ f%inodes2 = f%inodes2(:,new_index_list)
+
+
+ do e=1,f%nspec
+ do k=1,4
+ iglob = f%inodes1(k,e)
+ xyz(:,k) = nodes_coords(:,iglob)
+ enddo
+ xyz_c(1,e) = sum(xyz(1,:))/4
+ xyz_c(2,e) = sum(xyz(2,:))/4
+ xyz_c(3,e) = sum(xyz(3,:))/4
+ enddo
+
+
+ do e=1,f%nspec
+ do k=1,4
+ iglob = f%inodes2(k,e)
+ xyz(:,k) = nodes_coords(:,iglob)
+ enddo
+ xyz_c(1,e) = sum(xyz(1,:))/4
+ xyz_c(2,e) = sum(xyz(2,:))/4
+ xyz_c(3,e) = sum(xyz(3,:))/4
+ enddo
+
+
+ deallocate(xyz_c)
+
+ end subroutine reorder_fault_elements_single
+
+! ---------------------------------------------------------------------------------------------------
+ subroutine lex_order(xyz_c,loc,nnodes,nspec)
+
+ integer ,intent(in) :: nnodes,nspec
+ double precision, dimension(3,nspec) :: xyz_c
+
+ integer nglob !,nspec
+ integer iglob(nnodes),loc(nspec)
+ logical ifseg(nspec)
+ double precision xp(nspec),yp(nspec),zp(nspec)
+ double precision UTM_X_MIN,UTM_X_MAX
+
+ integer ispec,i,j,ier
+ integer ieoff,ilocnum,nseg,ioff,iseg,ig
+
+ integer, dimension(:), allocatable :: ind,ninseg,iwork
+ double precision, dimension(:), allocatable :: work
+
+ double precision :: SMALLVALTOL
+
+
+ xp=xyz_c(1,:)
+ yp=xyz_c(2,:)
+ zp=xyz_c(3,:)
+
+ UTM_X_MAX = maxval(xp)
+ UTM_X_MIN = minval(xp)
+ ! define geometrical tolerance based upon typical size of the model
+ SMALLVALTOL = 1.d-10 * dabs(UTM_X_MAX - UTM_X_MIN)
+
+ ! dynamically allocate arrays
+ allocate(ind(nspec), &
+ ninseg(nspec), &
+ iwork(nspec), &
+ work(nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays'
+
+
+ ! establish initial pointers
+ do ispec=1,nspec
+ ieoff=1*(ispec-1)
+ do ilocnum=1,1 !NGLLSQUARE
+ loc(ilocnum+ieoff)=ilocnum+ieoff
+ enddo
+ enddo
+
+ ifseg(:)=.false.
+
+ nseg=1
+ ifseg(1)=.true.
+ ninseg(1)=nspec
+
+ do j=1,NDIM
+
+ ! sort within each segment
+ ioff=1
+ do iseg=1,nseg
+ if(j == 1) then
+ call rank(xp(ioff),ind,ninseg(iseg))
+ else if(j == 2) then
+ call rank(yp(ioff),ind,ninseg(iseg))
+ else
+ call rank(zp(ioff),ind,ninseg(iseg))
+ endif
+ call swap_all(loc(ioff),xp(ioff),yp(ioff),zp(ioff),iwork,work,ind,ninseg(iseg))
+ ioff=ioff+ninseg(iseg)
+ enddo
+
+ ! check for jumps in current coordinate
+ ! compare the coordinates of the points within a small tolerance
+ if(j == 1) then
+ do i=2,nspec
+ if(dabs(xp(i)-xp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
+ enddo
+ else if(j == 2) then
+ do i=2,nspec
+ if(dabs(yp(i)-yp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
+ enddo
+ else
+ do i=2,nspec
+ if(dabs(zp(i)-zp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
+ enddo
+ endif
+
+ ! count up number of different segments
+ nseg=0
+ do i=1,nspec
+ if(ifseg(i)) then
+ nseg=nseg+1
+ ninseg(nseg)=1
+ else
+ ninseg(nseg)=ninseg(nseg)+1
+ endif
+ enddo
+ enddo
+
+ ! deallocate arrays
+ deallocate(ind)
+ deallocate(ninseg)
+ deallocate(iwork)
+ deallocate(work)
+
+ end subroutine lex_order
+
!===================================================================================================
!--------------------------------------------------
- ! Repartitioning : two coupled fault side1/side2 elements are transfered to the same partition
+ ! Repartitioning : two coupled faultside1/side2 elements are transfered to the same partition
!--------------------------------------------------
Subroutine fault_repartition (nelmnts, nnodes, elmnts, nsize, &
nproc, part, esize)
-! part(e) = index of the processor to which element #e is assigned.
-! Fault elements and neighbors are assigned to the same processor.
-! Part, once modified, will be input for write_partition_database.
+! part : iproc number of processor partioned. It will altered patching fault elements into the same partion.
+! Part, once is altered , will be input for write_partition_database.
!INPUTS
integer(long), intent(in) :: nelmnts,nsize
@@ -345,31 +499,49 @@
end subroutine fault_repartition
! ---------------------------------------------------------------------------------------------------
-! JPA parallel fault
-! subroutine fault_repartition_parallel (nelmnts, nnodes, elmnts, nsize, &
-! nproc, part, esize)
-!
-!!JPA loop over all faults
-!!JPA loop over all fault element pairs
-!!JPA assign the pair of elements to the processor with lowest rank among them
-!
-! do i=1,2 !JPA repeated to handle triple junctions
-!
-! do iflt=1,size(faults)
-! do e=1,faults(iflt)%nspec
-! e1 = faults(iflt)%ispec1(e)
-! e2 = faults(iflt)%ispec2(e)
-! proc1 = part(e1)
-! proc2 = part(e2)
-! part(e1) = min(proc1,proc2)
-! part(e2) = part(e1)
-! end do
-! end do
-!
-! end do
-!
-! end subroutine fault_repartition_parallel
+ subroutine fault_repartition_parallel (nelmnts, nnodes, elmnts, nsize, &
+ nproc, part, esize)
+
+
+!INPUTS
+ integer(long), intent(in) :: nelmnts,nsize
+ integer, intent(in) :: nnodes, nproc, esize
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+!OUTPUTS :
+ integer, dimension(0:nelmnts-1), intent(inout) :: part
+
+ integer, dimension(0:nnodes-1) :: nnodes_elmnts
+ integer, dimension(0:nsize*nnodes-1) :: nodes_elmnts
+ integer :: max_neighbour
+
+ integer :: i,j, ipart,nproc_null,nproc_null_final
+ integer :: k1, k2, k,e,iflt,inode
+ integer, dimension(:), allocatable :: elem_proc_null
+
+
+ integer :: e1,e2,proc1,proc2
+
+ do iflt=1,size(faults)
+ do e=1,faults(iflt)%nspec
+ e1 = faults(iflt)%ispec1(e) - 1
+ e2 = faults(iflt)%ispec2(e) - 1
+ proc1 = part(e1)
+ proc2 = part(e2)
+ part(e1) = min(proc1,proc2)
+ part(e2) = part(e1)
+ end do
+ end do
+
+ do iflt=1,size(faults)
+ do e=1,faults(iflt)%nspec
+ e1 = faults(iflt)%ispec1(e) - 1
+ e2 = faults(iflt)%ispec2(e) - 1
+ enddo
+ enddo
+
+ end subroutine fault_repartition_parallel
+
! ---------------------------------------------------------------------------------------------------
! See subroutine write_boundaries_database in part_decompose_mesh_SCOTCH.f90
!
@@ -400,6 +572,7 @@
! get number of fault elements in this partition
nspec_fault_1 = count( part(faults(iflt)%ispec1-1) == iproc )
nspec_fault_2 = count( part(faults(iflt)%ispec2-1) == iproc )
+
if (nspec_fault_1 /= nspec_fault_2) then
print *, 'Fault # ',iflt,', proc # ',iproc
print *, ' ispec1 : ', nspec_fault_1
@@ -448,4 +621,108 @@
end subroutine write_fault_database
+
+! -----------------------------------
+
+! sorting routines put in same file to allow for inlining
+
+ subroutine rank(A,IND,N)
+!
+! Use Heap Sort (Numerical Recipes)
+!
+ implicit none
+
+ integer n
+ double precision A(n)
+ integer IND(n)
+
+ integer i,j,l,ir,indx
+ double precision q
+
+ do j=1,n
+ IND(j)=j
+ enddo
+
+ if (n == 1) return
+
+ L=n/2+1
+ ir=n
+ 100 CONTINUE
+ IF (l>1) THEN
+ l=l-1
+ indx=ind(l)
+ q=a(indx)
+ ELSE
+ indx=ind(ir)
+ q=a(indx)
+ ind(ir)=ind(1)
+ ir=ir-1
+ if (ir == 1) then
+ ind(1)=indx
+ return
+ endif
+ ENDIF
+ i=l
+ j=l+l
+ 200 CONTINUE
+ IF (J <= IR) THEN
+ IF (J<IR) THEN
+ IF ( A(IND(j))<A(IND(j+1)) ) j=j+1
+ ENDIF
+ IF (q<A(IND(j))) THEN
+ IND(I)=IND(J)
+ I=J
+ J=J+J
+ ELSE
+ J=IR+1
+ ENDIF
+ goto 200
+ ENDIF
+ IND(I)=INDX
+ goto 100
+
+ end subroutine rank
+
+! ------------------------------------------------------------------
+
+ subroutine swap_all(IA,A,B,C,IW,W,ind,n)
+!
+! swap arrays IA, A, B and C according to addressing in array IND
+!
+ implicit none
+
+ integer n
+
+ integer IND(n)
+ integer IA(n),IW(n)
+ double precision A(n),B(n),C(n),W(n)
+
+ integer i
+
+ IW(:) = IA(:)
+ W(:) = A(:)
+
+ do i=1,n
+ IA(i)=IW(ind(i))
+ A(i)=W(ind(i))
+ enddo
+
+ W(:) = B(:)
+
+ do i=1,n
+ B(i)=W(ind(i))
+ enddo
+
+ W(:) = C(:)
+
+ do i=1,n
+ C(i)=W(ind(i))
+ enddo
+
+end subroutine swap_all
+
+! ------------------------------------------------------------------
+
+
end module fault_scotch
+
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90 2012-05-27 07:40:22 UTC (rev 20208)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90 2012-05-27 22:20:53 UTC (rev 20209)
@@ -56,7 +56,8 @@
use fault_object, only: fault_read_input,fault_setup, &
fault_save_arrays,fault_save_arrays_test,&
nnodes_coords_open,nodes_coords_open,ANY_FAULT_IN_THIS_PROC,&
- ANY_FAULT
+ ANY_FAULT
+ use fault_solve, only : PARALLEL_FAULT
implicit none
@@ -210,7 +211,7 @@
nodes_coords_ext_mesh,nnodes_ext_mesh,&
elmnts_ext_mesh,nelmnts_ext_mesh)
! at this point (xyz)store_dummy are still open
- call fault_setup (ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, &
+ if (.NOT.PARLLEL_FAULT) call fault_setup (ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, &
xstore,ystore,zstore,nspec,nglob,myrank)
! this closes (xyz)store_dummy
endif
@@ -229,6 +230,15 @@
num_interfaces_ext_mesh,max_interface_size_ext_mesh,&
my_neighbours_ext_mesh,NPROC)
+!SURENDRA (setting up parallel fault)
+ if (PARALLEL_FAULT .AND. ANY_FAULT) then
+ call sync_all()
+ !at this point (xyz)store_dummy are still open
+ call fault_setup (ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, &
+ xstore,ystore,zstore,nspec,nglob,myrank)
+ ! this closes (xyz)store_dummy
+ endif
+
! sets material velocities
call sync_all()
if( myrank == 0) then
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-27 07:40:22 UTC (rev 20208)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-05-27 22:20:53 UTC (rev 20209)
@@ -81,7 +81,7 @@
type bc_dynflt_type
private
- integer :: nspec,nglob
+ integer :: nspec,nglob,nglob_all,npoin_all
real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T0,T,V,D
real(kind=CUSTOM_REAL), dimension(:,:), pointer :: coord
real(kind=CUSTOM_REAL), dimension(:,:,:), pointer :: R
@@ -93,7 +93,8 @@
type(asperity_type), pointer :: asp => null()
logical :: allow_opening = .false. ! default : do not allow opening
type(dataT_type) :: dataT
- type(dataXZ_type) :: dataXZ
+ type(dataXZ_type) :: dataXZ,dataXZ_all
+ integer, dimension(:), pointer :: poin_offset,npoin_perproc
end type bc_dynflt_type
type(bc_dynflt_type), allocatable, save :: faults(:)
@@ -110,6 +111,8 @@
logical, save :: SIMULATION_TYPE_DYN = .false.
+ logical, save :: PARALLEL_FAULT = .false.
+
logical, save :: TPV16 = .false.
logical, save :: Rate_AND_State = .true.
@@ -129,12 +132,13 @@
! Minv inverse mass matrix
! dt global time step
!
-subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt)
+subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt,accel)
character(len=256), intent(in) :: prname ! 'proc***'
real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
double precision, intent(in) :: DTglobal
integer, intent(in) :: nt
+ real(kind=CUSTOM_REAL) :: accel(:,:)
real(kind=CUSTOM_REAL) :: dt
integer :: iflt,ier,dummy_idfault
@@ -187,7 +191,7 @@
allocate( faults(nbfaults) )
do iflt=1,nbfaults
read(IIN_PAR,nml=BEGIN_FAULT,end=100)
- call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
+ call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt,accel)
enddo
close(IIN_BIN)
close(IIN_PAR)
@@ -212,12 +216,14 @@
!---------------------------------------------------------------------
-subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
+subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt_tmp,NT,iflt,accel)
+ use specfem_par
+ real(kind=CUSTOM_REAL) :: accel(:,:)
type(bc_dynflt_type), intent(inout) :: bc
real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
integer, intent(in) :: IIN_BIN,IIN_PAR,NT,iflt
- real(kind=CUSTOM_REAL), intent(in) :: dt
+ real(kind=CUSTOM_REAL), intent(in) :: dt_tmp
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: jacobian2Dw
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: normal
@@ -270,23 +276,25 @@
NAMELIST / ASP / Fload,nFload
read(IIN_BIN) bc%nspec,bc%nglob
- if (bc%nspec==0) return
- allocate( bc%ibulk1(bc%nglob) )
- allocate( bc%ibulk2(bc%nglob) )
- allocate( ibool1(NGLLSQUARE,bc%nspec) )
- allocate(normal(NDIM,NGLLSQUARE,bc%nspec))
- allocate(jacobian2Dw(NGLLSQUARE,bc%nspec))
+ if (.NOT.PARALLEL_FAULT .and. bc%nspec==0) return
+ if (bc%nspec>0) then
+ allocate( bc%ibulk1(bc%nglob) )
+ allocate( bc%ibulk2(bc%nglob) )
+ allocate( ibool1(NGLLSQUARE,bc%nspec) )
+ allocate(normal(NDIM,NGLLSQUARE,bc%nspec))
+ allocate(jacobian2Dw(NGLLSQUARE,bc%nspec))
- allocate(bc%coord(3,(bc%nglob)))
- read(IIN_BIN) ibool1
- read(IIN_BIN) jacobian2Dw
- read(IIN_BIN) normal
- read(IIN_BIN) bc%ibulk1
- read(IIN_BIN) bc%ibulk2
- read(IIN_BIN) bc%coord(1,:)
- read(IIN_BIN) bc%coord(2,:)
- read(IIN_BIN) bc%coord(3,:)
- bc%dt = dt
+ allocate(bc%coord(3,(bc%nglob)))
+ read(IIN_BIN) ibool1
+ read(IIN_BIN) jacobian2Dw
+ read(IIN_BIN) normal
+ read(IIN_BIN) bc%ibulk1
+ read(IIN_BIN) bc%ibulk2
+ read(IIN_BIN) bc%coord(1,:)
+ read(IIN_BIN) bc%coord(2,:)
+ read(IIN_BIN) bc%coord(3,:)
+ endif
+ bc%dt = dt_tmp
allocate( bc%B(bc%nglob) )
bc%B = 0e0_CUSTOM_REAL
@@ -303,6 +311,63 @@
bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
enddo
enddo
+
+ if(PARALLEL_FAULT) then
+ accel=0._CUSTOM_REAL
+ accel(1,bc%ibulk1) = bc%B(:)
+ ! sends accel values to corresponding MPI interface neighbors
+ call assemble_MPI_vector_ext_mesh_s(NPROC,NGLOB_AB,accel, &
+ buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
+ my_neighbours_ext_mesh, &
+ request_send_vector_ext_mesh,request_recv_vector_ext_mesh)
+
+ ! assembles with other MPI processes
+ call assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,accel, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ my_neighbours_ext_mesh)
+
+ ! waits for send/receive requests to be completed and assembles values
+ call assemble_MPI_vector_ext_mesh_w(NPROC,NGLOB_AB,accel, &
+ buffer_recv_vector_ext_mesh,num_interfaces_ext_mesh,&
+ max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ request_send_vector_ext_mesh,request_recv_vector_ext_mesh)
+
+ bc%B(:) = accel(1,bc%ibulk1)
+
+ accel=0._CUSTOM_REAL
+ accel(1,bc%ibulk1) = nx(:)
+ accel(2,bc%ibulk1) = ny(:)
+ accel(3,bc%ibulk1) = nz(:)
+ ! sends accel values to corresponding MPI interface neighbors
+ call assemble_MPI_vector_ext_mesh_s(NPROC,NGLOB_AB,accel, &
+ buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
+ my_neighbours_ext_mesh, &
+ request_send_vector_ext_mesh,request_recv_vector_ext_mesh)
+
+ ! assembles with other MPI processes
+ call assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,accel, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ my_neighbours_ext_mesh)
+
+ ! waits for send/receive requests to be completed and assembles values
+ call assemble_MPI_vector_ext_mesh_w(NPROC,NGLOB_AB,accel, &
+ buffer_recv_vector_ext_mesh,num_interfaces_ext_mesh,&
+ max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ request_send_vector_ext_mesh,request_recv_vector_ext_mesh)
+
+ nx(:) = accel(1,bc%ibulk1)
+ ny(:) = accel(2,bc%ibulk1)
+ nz(:) = accel(3,bc%ibulk1)
+ endif
+
do k=1,bc%nglob
norm = sqrt( nx(k)*nx(k) + ny(k)*ny(k) + nz(k)*nz(k) )
nx(k) = nx(k) / norm
@@ -711,7 +776,7 @@
if (.not. allocated(faults)) return
do iflt=1,size(faults)
- if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d(faults(iflt),F,Vel,Dis,iflt)
+ call BC_DYNFLT_set3d(faults(iflt),F,Vel,Dis,iflt)
enddo
end subroutine bc_dynflt_set3d_all
@@ -719,7 +784,7 @@
!---------------------------------------------------------------------
subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
- use specfem_par, only:it,NSTEP
+ use specfem_par, only:it,NSTEP,myrank
real(kind=CUSTOM_REAL), intent(inout) :: MxA(:,:)
type(bc_dynflt_type), intent(inout) :: bc
@@ -740,155 +805,166 @@
real(kind=CUSTOM_REAL) :: time
real(kind=CUSTOM_REAL) :: psi
- half_dt = 0.5e0_CUSTOM_REAL*bc%dt
- Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+ if(bc%nspec > 0) then !Surendra : for parallel faults
- ! get predicted values
- dD = get_jump(bc,D) ! dD_predictor
- dV = get_jump(bc,V) ! dV_predictor
- dA = get_weighted_jump(bc,MxA) ! dA_free
+ half_dt = 0.5e0_CUSTOM_REAL*bc%dt
+ Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+
+ ! get predicted values
+ dD = get_jump(bc,D) ! dD_predictor
+ dV = get_jump(bc,V) ! dV_predictor
+ dA = get_weighted_jump(bc,MxA) ! dA_free
- ! rotate to fault frame (tangent,normal)
- ! component 3 is normal to the fault
- dD = rotate(bc,dD,1)
- dV = rotate(bc,dV,1)
- dA = rotate(bc,dA,1)
+ ! rotate to fault frame (tangent,normal)
+ ! component 3 is normal to the fault
+ dD = rotate(bc,dD,1)
+ dV = rotate(bc,dV,1)
+ dA = rotate(bc,dA,1)
+
+ ! T_stick
+ T(1,:) = bc%Z * ( dV(1,:) + half_dt*dA(1,:) )
+ T(2,:) = bc%Z * ( dV(2,:) + half_dt*dA(2,:) )
+ T(3,:) = bc%Z * ( dV(3,:) + half_dt*dA(3,:) )
+
+ !Warning : dirty particular free surface condition z = 0.
+ ! where (bc%zcoord(:) > - SMALLVAL) T(2,:) = 0
+ ! do k=1,bc%nglob
+ ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) < SMALLVAL) T(2,k) = 0.e0_CUSTOM_REAL
+ ! end do
- ! T_stick
- T(1,:) = bc%Z * ( dV(1,:) + half_dt*dA(1,:) )
- T(2,:) = bc%Z * ( dV(2,:) + half_dt*dA(2,:) )
- T(3,:) = bc%Z * ( dV(3,:) + half_dt*dA(3,:) )
+ ! add initial stress
+ T = T + bc%T0
+
+ ! Solve for normal stress (negative is compressive)
+ ! Opening implies free stress
+ if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
+
+ if(.not. Rate_AND_State) then ! Update slip weakening friction:
+ ! Update slip state variable
+ ! WARNING: during opening the friction state variable should not evolve
+ theta_old = bc%swf%theta
+ call swf_update_state(bc%D,dD,bc%V,bc%swf)
+
+ ! Update friction coeficient
+ bc%MU = swf_mu(bc%swf)
+
+ ! combined with time-weakening for nucleation
+ ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
+ if (TPV16) then
+ where (bc%swf%T <= it*bc%dt) bc%MU = bc%swf%mud
+ endif
- !Warning : dirty particular free surface condition z = 0.
- ! where (bc%zcoord(:) > - SMALLVAL) T(2,:) = 0
- ! do k=1,bc%nglob
- ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) < SMALLVAL) T(2,k) = 0.e0_CUSTOM_REAL
- ! end do
-
- ! add initial stress
- T = T + bc%T0
-
- ! Solve for normal stress (negative is compressive)
- ! Opening implies free stress
- if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
-
- if(.not. Rate_AND_State) then ! Update slip weakening friction:
- ! Update slip state variable
- ! WARNING: during opening the friction state variable should not evolve
- theta_old = bc%swf%theta
- call swf_update_state(bc%D,dD,bc%V,bc%swf)
-
- ! Update friction coeficient
- bc%MU = swf_mu(bc%swf)
-
- ! combined with time-weakening for nucleation
- ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
- if (TPV16) then
- where (bc%swf%T <= it*bc%dt) bc%MU = bc%swf%mud
+ ! Update strength
+ strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
+
+ ! Solve for shear stress
+ tStick = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
+ tnew = min(tStick,strength)
+ tStick = max(tStick,1e0_CUSTOM_REAL)
+ T(1,:) = tnew * T(1,:)/tStick
+ T(2,:) = tnew * T(2,:)/tStick
+
+ else ! Update rate and state friction:
+
+ ! smooth loading within nucleation patch
+ !WARNING : ad hoc for SCEC benchmark TPV10x
+ TxExt = 0._CUSTOM_REAL
+ TLoad = 1.0_CUSTOM_REAL
+ DTau0 = 25e6_CUSTOM_REAL
+ time = it*bc%dt !time will never be zero. it starts from 1
+ if (time <= TLoad) then
+ GLoad = exp( (time-TLoad)*(time-Tload) / (time*(time-2.0_CUSTOM_REAL*TLoad)) )
+ else
+ GLoad = 1.0_CUSTOM_REAL
+ endif
+ TxExt = DTau0 * bc%asp%Fload * GLoad
+ T(1,:) = T(1,:) + TxExt
+
+ !JPA the solver below can be refactored into a loop with two passes
+ Vf = Vnorm_old
+ theta_old = bc%rsf%theta
+ call rsf_update_state(Vf,bc%dt,bc%rsf)
+
+ tStick = sqrt(T(1,:)**2 + T(2,:)**2)
+ do iNode=1,bc%nglob
+ Vf1(iNode)=rtsafe(funcd,0.0,Vf(iNode)+5.0,1e-5,tStick(iNode),-T(3,iNode),bc%Z(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),bc%rsf%L(iNode),bc%rsf%theta(iNode))
+ tau1(iNode) = tStick(iNode) - bc%Z(iNode)*Vf1(iNode)
+ enddo
+
+ ! Updating state variable: 2nd loop
+ Vfavg = 0.5e0_CUSTOM_REAL*(Vf + Vf1)
+ bc%rsf%theta = theta_old
+ call rsf_update_state(Vfavg,bc%dt,bc%rsf)
+
+ ! NR search 2nd loop
+ do iNode=1,bc%nglob
+ Vf2(iNode)=rtsafe(funcd,0.0,Vf(iNode)+5.0,1e-5,tStick(iNode),-T(3,iNode),bc%Z(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),bc%rsf%L(iNode),bc%rsf%theta(iNode))
+ tau2(iNode) = tStick(iNode) - bc%Z(iNode)*Vf2(iNode)
+ enddo
+
+ tStick = max(tStick,1e0_CUSTOM_REAL)
+ T(1,:) = tau2 * T(1,:)/tStick
+ T(2,:) = tau2 * T(2,:)/tStick
+
endif
+
+ ! Save total tractions
+ bc%T = T
+
+ ! Subtract initial stress
+ T = T - bc%T0
+
+ ! Update slip acceleration da=da_free-T/(0.5*dt*Z)
+ dA(1,:) = dA(1,:) - T(1,:)/(bc%Z*half_dt)
+ dA(2,:) = dA(2,:) - T(2,:)/(bc%Z*half_dt)
+ dA(3,:) = dA(3,:) - T(3,:)/(bc%Z*half_dt)
+
+ ! Update slip and slip rate, in fault frame
+ bc%D = dD
+ bc%V = dV + half_dt*dA
+
+ ! Rotate tractions back to (x,y,z) frame
+ T = rotate(bc,T,-1)
+
+ ! Add boundary term B*T to M*a
+ MxA(1,bc%ibulk1) = MxA(1,bc%ibulk1) + bc%B*T(1,:)
+ MxA(2,bc%ibulk1) = MxA(2,bc%ibulk1) + bc%B*T(2,:)
+ MxA(3,bc%ibulk1) = MxA(3,bc%ibulk1) + bc%B*T(3,:)
- ! Update strength
- strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
+ MxA(1,bc%ibulk2) = MxA(1,bc%ibulk2) - bc%B*T(1,:)
+ MxA(2,bc%ibulk2) = MxA(2,bc%ibulk2) - bc%B*T(2,:)
+ MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
- ! Solve for shear stress
- tStick = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
- tnew = min(tStick,strength)
- tStick = max(tStick,1e0_CUSTOM_REAL)
- T(1,:) = tnew * T(1,:)/tStick
- T(2,:) = tnew * T(2,:)/tStick
-
- else ! Update rate and state friction:
-
- ! smooth loading within nucleation patch
- !WARNING : ad hoc for SCEC benchmark TPV10x
- TxExt = 0._CUSTOM_REAL
- TLoad = 1.0_CUSTOM_REAL
- DTau0 = 25e6_CUSTOM_REAL
- time = it*bc%dt !time will never be zero. it starts from 1
- if (time <= TLoad) then
- GLoad = exp( (time-TLoad)*(time-Tload) / (time*(time-2.0_CUSTOM_REAL*TLoad)) )
+ !-- intermediate storage of outputs --
+ Vnorm = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+ if(.not. Rate_AND_State) then
+ theta_new = bc%swf%theta
+ dc = bc%swf%dc
else
- GLoad = 1.0_CUSTOM_REAL
+ theta_new = bc%rsf%theta
+ dc = bc%rsf%L
endif
- TxExt = DTau0 * bc%asp%Fload * GLoad
- T(1,:) = T(1,:) + TxExt
+
+ call store_dataXZ(bc%dataXZ, strength, theta_old, theta_new, dc, &
+ Vnorm_old, Vnorm, it*bc%dt,bc%dt)
+ call store_dataT(bc%dataT,bc%D,bc%V,bc%T,theta_new,it)
- !JPA the solver below can be refactored into a loop with two passes
- Vf = Vnorm_old
- theta_old = bc%rsf%theta
- call rsf_update_state(Vf,bc%dt,bc%rsf)
-
- tStick = sqrt(T(1,:)**2 + T(2,:)**2)
- do iNode=1,bc%nglob
- Vf1(iNode)=rtsafe(funcd,0.0,Vf(iNode)+5.0,1e-5,tStick(iNode),-T(3,iNode),bc%Z(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),bc%rsf%L(iNode),bc%rsf%theta(iNode))
- tau1(iNode) = tStick(iNode) - bc%Z(iNode)*Vf1(iNode)
- enddo
-
- ! Updating state variable: 2nd loop
- Vfavg = 0.5e0_CUSTOM_REAL*(Vf + Vf1)
- bc%rsf%theta = theta_old
- call rsf_update_state(Vfavg,bc%dt,bc%rsf)
-
- ! NR search 2nd loop
- do iNode=1,bc%nglob
- Vf2(iNode)=rtsafe(funcd,0.0,Vf(iNode)+5.0,1e-5,tStick(iNode),-T(3,iNode),bc%Z(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),bc%rsf%L(iNode),bc%rsf%theta(iNode))
- tau2(iNode) = tStick(iNode) - bc%Z(iNode)*Vf2(iNode)
- enddo
-
- tStick = max(tStick,1e0_CUSTOM_REAL)
- T(1,:) = tau2 * T(1,:)/tStick
- T(2,:) = tau2 * T(2,:)/tStick
-
+ !-- outputs --
+ ! write dataT every NTOUT time step or at the end of simulation
+ if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it)
+ if ( it == NSTEP) call SCEC_Write_RuptureTime(bc%dataXZ,bc%dt,NSTEP,iflt)
endif
- ! Save total tractions
- bc%T = T
-
- ! Subtract initial stress
- T = T - bc%T0
-
- ! Update slip acceleration da=da_free-T/(0.5*dt*Z)
- dA(1,:) = dA(1,:) - T(1,:)/(bc%Z*half_dt)
- dA(2,:) = dA(2,:) - T(2,:)/(bc%Z*half_dt)
- dA(3,:) = dA(3,:) - T(3,:)/(bc%Z*half_dt)
-
- ! Update slip and slip rate, in fault frame
- bc%D = dD
- bc%V = dV + half_dt*dA
-
- ! Rotate tractions back to (x,y,z) frame
- T = rotate(bc,T,-1)
-
- ! Add boundary term B*T to M*a
- MxA(1,bc%ibulk1) = MxA(1,bc%ibulk1) + bc%B*T(1,:)
- MxA(2,bc%ibulk1) = MxA(2,bc%ibulk1) + bc%B*T(2,:)
- MxA(3,bc%ibulk1) = MxA(3,bc%ibulk1) + bc%B*T(3,:)
-
- MxA(1,bc%ibulk2) = MxA(1,bc%ibulk2) - bc%B*T(1,:)
- MxA(2,bc%ibulk2) = MxA(2,bc%ibulk2) - bc%B*T(2,:)
- MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
-
- !-- intermediate storage of outputs --
- Vnorm = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
- if(.not. Rate_AND_State) then
- theta_new = bc%swf%theta
- dc = bc%swf%dc
- else
- theta_new = bc%rsf%theta
- dc = bc%rsf%L
+ ! write dataXZ every NSNAP time step
+ if ( mod(it,NSNAP) == 0) then
+ if(PARALLEL_FAULT) then
+ call write_dataXZ(bc%dataXZ,it,iflt)
+ else
+ call gather_dataXZ(bc)
+ if(myrank==0) call write_dataXZ(bc%dataXZ_all,it,iflt)
+ endif
endif
- call store_dataXZ(bc%dataXZ, strength, theta_old, theta_new, dc, &
- Vnorm_old, Vnorm, it*bc%dt,bc%dt)
- call store_dataT(bc%dataT,bc%D,bc%V,bc%T,theta_new,it)
-
- !-- outputs --
- ! write dataT every NTOUT time step or at the end of simulation
- if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it)
- ! write dataXZ every NSNAP time step
- if ( mod(it,NSNAP) == 0) call write_dataXZ(bc%dataXZ,it,iflt)
- if ( it == NSTEP) call SCEC_Write_RuptureTime(bc%dataXZ,bc%dt,NSTEP,iflt)
-
end subroutine BC_DYNFLT_set3d
!===============================================================
@@ -1315,10 +1391,17 @@
!-------------------------------------------------------------------------------------------------
subroutine init_dataXZ(DataXZ,bc,nglob)
+ use specfem_par, only : NPROC,myrank
+
type(dataXZ_type), intent(inout) :: DataXZ
type(bc_dynflt_type) :: bc
integer, intent(in) :: nglob
+ integer :: npoin_all,npoin,iproc,ier
+
+ DataXZ%npoin = nglob !Surendra
+ npoin = DataXZ%npoin
+
allocate(DataXZ%stg(nglob))
if(.not. Rate_AND_State) then
DataXZ%sta => bc%swf%theta
@@ -1342,8 +1425,66 @@
DataXZ%tRUP = 0e0_CUSTOM_REAL
DataXZ%tPZ = 0e0_CUSTOM_REAL
+ !Surendra : for parallel fault
+ if(PARLLEL_FAULT) then
+ call sum_all_i(bc%nglob,bc%nglob_all)
+ if(myrank==0) then
+ npoin_all=bc%nglob_all
+ !bc%DataXZ_all%npoin=npoin_all
+ bc%npoin_all=npoin_all
+ allocate(bc%DataXZ_all%xcoord(npoin_all))
+ allocate(bc%DataXZ_all%ycoord(npoin_all))
+ allocate(bc%DataXZ_all%zcoord(npoin_all))
+ allocate(bc%DataXZ_all%t1(npoin_all))
+ allocate(bc%DataXZ_all%t2(npoin_all))
+ allocate(bc%DataXZ_all%t3(npoin_all))
+ allocate(bc%DataXZ_all%d1(npoin_all))
+ allocate(bc%DataXZ_all%d2(npoin_all))
+ allocate(bc%DataXZ_all%v1(npoin_all))
+ allocate(bc%DataXZ_all%v2(npoin_all))
+ allocate(bc%DataXZ_all%tRUP(npoin_all))
+ allocate(bc%DataXZ_all%tPZ(npoin_all))
+ allocate(bc%DataXZ_all%stg(npoin_all))
+ allocate(bc%DataXZ_all%sta(npoin_all))
+ endif
+
+ allocate(bc%npoin_perproc(NPROC))
+ allocate(bc%poin_offset(NPROC))
+ bc%npoin_perproc=0
+ call gather_all_i(npoin,1,bc%npoin_perproc,1,NPROC)
+
+ bc%poin_offset(1)=0
+ do iproc=2,NPROC
+ bc%poin_offset(iproc) = sum(bc%npoin_perproc(1:iproc-1))
+ enddo
+
+ call gatherv_all_cr(DataXZ%xcoord,DataXZ%npoin,bc%DataXZ_all%xcoord,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ call gatherv_all_cr(DataXZ%ycoord,DataXZ%npoin,bc%DataXZ_all%ycoord,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ call gatherv_all_cr(DataXZ%zcoord,DataXZ%npoin,bc%DataXZ_all%zcoord,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ endif
+
end subroutine init_dataXZ
+!---------------------------------------------------------------
+subroutine gather_dataXZ(bc)
+ use specfem_par, only : NPROC
+
+ type(bc_dynflt_type) :: bc
+
+ call gatherv_all_cr(bc%DataXZ%t1,bc%DataXZ%npoin,bc%DataXZ_all%t1,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ call gatherv_all_cr(bc%DataXZ%t2,bc%DataXZ%npoin,bc%DataXZ_all%t2,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ call gatherv_all_cr(bc%DataXZ%t3,bc%DataXZ%npoin,bc%DataXZ_all%t3,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ call gatherv_all_cr(bc%DataXZ%d1,bc%DataXZ%npoin,bc%DataXZ_all%d1,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ call gatherv_all_cr(bc%DataXZ%d2,bc%DataXZ%npoin,bc%DataXZ_all%d2,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ call gatherv_all_cr(bc%DataXZ%v1,bc%DataXZ%npoin,bc%DataXZ_all%v1,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ call gatherv_all_cr(bc%DataXZ%v2,bc%DataXZ%npoin,bc%DataXZ_all%v2,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ call gatherv_all_cr(bc%DataXZ%tRUP,bc%DataXZ%npoin,bc%DataXZ_all%tRUP,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ call gatherv_all_cr(bc%DataXZ%tPZ,bc%DataXZ%npoin,bc%DataXZ_all%tPZ,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ call gatherv_all_cr(bc%DataXZ%stg,bc%DataXZ%npoin,bc%DataXZ_all%stg,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+ call gatherv_all_cr(bc%DataXZ%sta,bc%DataXZ%npoin,bc%DataXZ_all%sta,bc%npoin_perproc,bc%poin_offset,bc%npoin_all,NPROC)
+
+end subroutine gather_dataXZ
+
!---------------------------------------------------------------
subroutine store_dataXZ(dataXZ,stg,dold,dnew,dc,vold,vnew,time,dt)
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/iterate_time.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/iterate_time.f90 2012-05-27 07:40:22 UTC (rev 20208)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/iterate_time.f90 2012-05-27 22:20:53 UTC (rev 20209)
@@ -64,7 +64,6 @@
do it = 1,NSTEP
- write(IMAIN,*) 'it:',it !test
! simulation status output and stability check
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5) then
call it_check_stability()
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/prepare_timerun.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/prepare_timerun.f90 2012-05-27 07:40:22 UTC (rev 20208)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/prepare_timerun.f90 2012-05-27 22:20:53 UTC (rev 20209)
@@ -115,7 +115,7 @@
call prepare_timerun_mass_matrices()
! Loading kinematic and dynamic fault solvers.
- call BC_DYNFLT_init(prname,rmass,DT,NSTEP)
+ call BC_DYNFLT_init(prname,rmass,DT,NSTEP,accel)
call BC_KINFLT_init(prname,rmass,DT,NSTEP)
More information about the CIG-COMMITS
mailing list