Skip to content

Commit

Permalink
write out moment rate function as stf
Browse files Browse the repository at this point in the history 
verified in homogeneous ridgecrest model
  • Loading branch information
@elifo
elifo committed Jul 1, 2022
1 parent 0e1a590 commit 77abef3
Show file tree
Hide file tree
Showing 4 changed files with 114 additions and 39 deletions.
5 changes: 4 additions & 1 deletion src/generate_databases/create_regions_mesh.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ subroutine create_regions_mesh()
fault_save_arrays,fault_save_arrays_test, &
fault_save_arrays_txt, &
nnodes_coords_open,nodes_coords_open,ANY_FAULT_IN_THIS_PROC, &
ANY_FAULT
ANY_FAULT, fault_set_material

implicit none

Expand Down Expand Up @@ -204,6 +204,9 @@ subroutine create_regions_mesh()
endif
call get_model()

! store fault mu
call fault_set_material(nspec,mustore)

! sets up acoustic-elastic-poroelastic coupling surfaces
call synchronize_all()
if (myrank == 0) then
Expand Down
78 changes: 73 additions & 5 deletions src/generate_databases/fault_generate_databases.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,8 @@ module fault_generate_databases
type fault_db_type
private
real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoordbulk1,ycoordbulk1,zcoordbulk1, &
xcoordbulk2,ycoordbulk2,zcoordbulk2
xcoordbulk2,ycoordbulk2,zcoordbulk2, &
mubulk1, mubulk2
real(kind=CUSTOM_REAL), dimension(:,:), pointer:: jacobian2Dw
real(kind=CUSTOM_REAL), dimension(:,:,:), pointer:: normal
real(kind=CUSTOM_REAL) :: eta
Expand Down Expand Up @@ -85,7 +86,7 @@ module fault_generate_databases
iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces

public :: fault_read_input, fault_setup, fault_save_arrays_test, fault_save_arrays, fault_save_arrays_txt, &
nodes_coords_open, nnodes_coords_open, ANY_FAULT_IN_THIS_PROC, ANY_FAULT
nodes_coords_open, nnodes_coords_open, ANY_FAULT_IN_THIS_PROC, ANY_FAULT, fault_set_material

contains

Expand Down Expand Up @@ -520,13 +521,13 @@ end subroutine close_fault

subroutine save_fault_xyzcoord_ibulk(fdb)

use create_regions_mesh_ext_par, only: xstore_unique,ystore_unique,zstore_unique
use create_regions_mesh_ext_par, only: xstore_unique,ystore_unique,zstore_unique,mustore

implicit none
type(fault_db_type), intent(inout) :: fdb

! local parameters
integer :: K1, K2, i, ier
integer :: K1, K2, i, ier, e, ie, je, ke, k

allocate( fdb%xcoordbulk1(fdb%nglob) ,stat=ier)
if (ier /= 0) call exit_MPI_without_rank('error allocating array 883')
Expand All @@ -540,6 +541,10 @@ subroutine save_fault_xyzcoord_ibulk(fdb)
if (ier /= 0) call exit_MPI_without_rank('error allocating array 887')
allocate( fdb%zcoordbulk2(fdb%nglob) ,stat=ier)
if (ier /= 0) call exit_MPI_without_rank('error allocating array 888')
allocate( fdb%mubulk1(fdb%nglob) ,stat=ier)
if (ier /= 0) call exit_MPI_without_rank('error allocating array 889')
allocate( fdb%mubulk2(fdb%nglob) ,stat=ier)
if (ier /= 0) call exit_MPI_without_rank('error allocating array 890')
fdb%xcoordbulk1(:) = 0.0_CUSTOM_REAL; fdb%ycoordbulk1(:) = 0.0_CUSTOM_REAL; fdb%zcoordbulk1(:) = 0.0_CUSTOM_REAL
fdb%xcoordbulk2(:) = 0.0_CUSTOM_REAL; fdb%ycoordbulk2(:) = 0.0_CUSTOM_REAL; fdb%zcoordbulk2(:) = 0.0_CUSTOM_REAL

Expand All @@ -555,6 +560,8 @@ subroutine save_fault_xyzcoord_ibulk(fdb)
fdb%zcoordbulk2(i) = zstore_unique(K2)
enddo

enddo

end subroutine save_fault_xyzcoord_ibulk


Expand Down Expand Up @@ -722,6 +729,7 @@ subroutine save_one_fault_test(f,IOUT)
write(IOUT,*) f%ibulk2(k),f%xcoordbulk2(k),f%ycoordbulk2(k),f%zcoordbulk2(k)
enddo


end subroutine save_one_fault_test

!=================================================================================
Expand Down Expand Up @@ -806,6 +814,10 @@ subroutine save_one_fault_bin(f,IOUT)
write(IOUT) f%ispec1
write(IOUT) f%ispec2

! Elif: store mu values
write(IOUT) f%mubulk1
write(IOUT) f%mubulk2

end subroutine save_one_fault_bin

!====================================================================================
Expand Down Expand Up @@ -864,10 +876,66 @@ subroutine save_one_fault_test_txt(f,IOUT)
write(IOUT,*) f%ibulk2(k),f%xcoordbulk2(k),f%ycoordbulk2(k),f%zcoordbulk2(k)
enddo

end subroutine save_one_fault_test_txt
end subroutine save_one_fault_test_txt
!====================================================================================

subroutine fault_set_material(nspec,mustore)

implicit none

integer, intent(in) :: nspec
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(in) :: mustore

! local parameters
integer :: iflt

! only partitions with a fault need to setup
if (.not. ANY_FAULT_IN_THIS_PROC) return

do iflt = 1,size(fault_db)

! checks if anything to do on this fault
if (fault_db(iflt)%nspec == 0) cycle

call save_fault_mu(fault_db(iflt), nspec, mustore)

enddo

end subroutine fault_set_material
!====================================================================================

subroutine save_fault_mu(fdb, nspec, mustore)

implicit none

type(fault_db_type), intent(inout) :: fdb
integer, intent(in) :: nspec
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(in) :: mustore

integer :: e, k, ie, je, ke, K1, K2

write(*,*) 'ELIF: fault_set_mu: ', size(mustore,1)
do e = 1, fdb%nspec
do k = 1, NGLLSQUARE
ie = fdb%ijk1(1,k,e)
je = fdb%ijk1(2,k,e)
ke = fdb%ijk1(3,k,e)
K1 = fdb%ibool1(k,e)
fdb%mubulk1(K1) = mustore(ie,je,ke,fdb%ispec1(e))

ie = fdb%ijk2(1,k,e)
je = fdb%ijk2(2,k,e)
ke = fdb%ijk2(3,k,e)
K2 = fdb%ibool2(k,e)
fdb%mubulk2(K2) = mustore(ie,je,ke,fdb%ispec2(e))

enddo
enddo

! write(*,*) 'ELIF: mubulk1: ', fdb%mubulk1(1:50)
! write(*,*) 'ELIF: mubulk2: ', fdb%mubulk1(1:50)

end subroutine save_fault_mu
!====================================================================================

end module fault_generate_databases
52 changes: 32 additions & 20 deletions src/specfem3D/fault_solver_common.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ module fault_solver_common
real(kind=CUSTOM_REAL) :: dt
integer, dimension(:), pointer :: ibulk1 => null(), ibulk2 => null() ! global nodes id
integer, dimension(:), pointer :: ispec1 => null(), ispec2 => null() ! fault elements id
real(kind=CUSTOM_REAL), dimension(:), pointer :: mubulk1 => null(), mubulk2 => null() ! fault mu
end type fault_type

! outputs(dyn) /inputs (kind) at selected times for all fault nodes:
Expand Down Expand Up @@ -205,6 +206,11 @@ subroutine initialize_fault (bc,IIN_BIN)
! fault elements ispec (touching lower surface)
allocate(bc%ispec2(bc%nspec),stat=ier)
if (ier /= 0) call exit_MPI_without_rank('error allocating array 2171')
! fault rigidity mu=rho*Vs*Vs
allocate(bc%mubulk1(bc%nglob),stat=ier)
if (ier /= 0) call exit_MPI_without_rank('error allocating array 2172')
allocate(bc%mubulk2(bc%nglob),stat=ier)
if (ier /= 0) call exit_MPI_without_rank('error allocating array 2173')

read(IIN_BIN) ibool1
read(IIN_BIN) jacobian2Dw
Expand All @@ -225,6 +231,10 @@ subroutine initialize_fault (bc,IIN_BIN)
! dummy value if not stored
if (ier /= 0) bc%ispec2(:) = 0

! Elif: read mu stored
read(IIN_BIN) bc%mubulk1
read(IIN_BIN) bc%mubulk2

! done reading faults_db.bin file, no more records in read(IIN_BIN).. from here on

! sets simulation time step size
Expand Down Expand Up @@ -774,41 +784,43 @@ subroutine SCEC_write_dataT(dataT)
end subroutine SCEC_write_dataT

!------------------------------------------------------------------------

subroutine write_STF_GPU(STF, NT, NGLOB)
! Write out source time function (Elif, 06/22)
! I could not assemble STF matrix of all procs
! Current version writes one file per proc.
subroutine write_STF_GPU(B, mu1, mu2, STF, NT, NGLOB, myrank)

use specfem_par, only: OUTPUT_FILES

implicit none

integer, intent(in) :: NT, NGLOB
integer, intent(in) :: NT, NGLOB, myrank
real(kind=CUSTOM_REAL), dimension(NGLOB,NT),intent(in) :: STF
real(kind=CUSTOM_REAL), dimension(NGLOB),intent(in) :: B, mu1, mu2
character(len=5):: iproc

! local parameters
integer :: i,iol
!real(kind=CUSTOM_REAL) :: moment_rate(NT)
integer :: i
integer, parameter :: IOUT_SC = 121 !WARNING: not very robust.
real(kind=CUSTOM_REAL) :: mu(NGLOB)

write(*,*) 'ELIF::write_STF_GPU, NT,NGLOB: ', NT, NGLOB
write(*,*) size(STF,1), size(STF,2)
write(*,*) 'ELIF::write_STF_GPU, NT,NGLOB, myrank: ', NT, NGLOB, myrank
write(*,*) 'size(STF,1), size(STF,2): ', size(STF,1), size(STF,2)
write(*,*) 'size(B,1) : ', size(B,1)
write(*,*) 'size(mu,1) : ', size(mu1,1)
write(*,*) 'maxval(mu1), maxval(mu2): ', maxval(mu1),maxval(mu2)

! write out
! open(IOUT_SC,file=trim(OUTPUT_FILES)//'STF.dat',status='replace')
! do i=1,NT
! write(IOUT_SC,*) STF(:,i)
! enddo
! close(IOUT_SC)

! bin file
INQUIRE( IOLENGTH=iol ) STF(:,1)
write(*,*) 'iol: ', iol
open(IOUT_SC,file=trim(OUTPUT_FILES)//'STF.dat',status='replace',access='direct',recl=iol)
! to verify what value to use when aysmmetrical fault
! take ave. for now
mu = (mu1+ mu2)* 0.5_CUSTOM_REAL

! STF = Moment rate (time) = sum(mu*B*Vslip)
write(iproc,'(i5.5)') myrank
open(IOUT_SC,file=trim(OUTPUT_FILES)//'proc_'//trim(iproc)//'_STF.dat',status='unknown')
do i=1,NT
write(IOUT_SC,rec=i) STF(:,i)
write(IOUT_SC,*) sum( mu(:)* B(:)* STF(:,i) )
enddo
close(IOUT_SC)


end subroutine write_STF_GPU
!------------------------------------------------------------------------

Expand Down
18 changes: 5 additions & 13 deletions src/specfem3D/fault_solver_dynamic.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2051,6 +2051,7 @@ subroutine init_dataXZ(dataXZ,bc)

if ( GPU_MODE) then
dataXZ%tRUP => bc%Trup(1,:)
dataXZ%STF => bc%STF
else
allocate(dataXZ%tRUP(bc%nglob),stat=ier)
if (ier /= 0) call exit_MPI_without_rank('error allocating array 1402')
Expand Down Expand Up @@ -2139,7 +2140,8 @@ subroutine init_dataXZ(dataXZ,bc)
bc%dataXZ_all%tRUP(1), &
bc%dataXZ_all%tPZ(1), &
bc%dataXZ_all%stg(1), &
bc%dataXZ_all%sta(1))
bc%dataXZ_all%sta(1),&
bc%dataXZ_all%STF(1,1))
endif

!note: crayftn compiler warns about possible copy which may slow down the code for dataXZ%npoin,dataXZ%xcoord,..
Expand Down Expand Up @@ -2180,7 +2182,7 @@ subroutine gather_dataXZ(bc)
implicit none

type(bc_dynandkinflt_type), intent(inout) :: bc
integer :: i
!integer :: i

! collects data from all processes onto main process arrays
call gatherv_all_cr(bc%dataXZ%t1,bc%dataXZ%npoin,bc%dataXZ_all%t1,bc%npoin_perproc,bc%poin_offset,bc%dataXZ_all%npoin,NPROC)
Expand All @@ -2194,16 +2196,6 @@ subroutine gather_dataXZ(bc)
call gatherv_all_cr(bc%dataXZ%tPZ,bc%dataXZ%npoin,bc%dataXZ_all%tPZ,bc%npoin_perproc,bc%poin_offset,bc%dataXZ_all%npoin,NPROC)
call gatherv_all_cr(bc%dataXZ%stg,bc%dataXZ%npoin,bc%dataXZ_all%stg,bc%npoin_perproc,bc%poin_offset,bc%dataXZ_all%npoin,NPROC)
call gatherv_all_cr(bc%dataXZ%sta,bc%dataXZ%npoin,bc%dataXZ_all%sta,bc%npoin_perproc,bc%poin_offset,bc%dataXZ_all%npoin,NPROC)
if ( size(bc%STF,1) > 0 ) then
write(*,*) 'size(bc%STF), bc%dataXZ_all%STF', size(bc%STF,1), size(bc%STF,2), &
size(bc%dataXZ_all%STF,1) , &
size(bc%dataXZ_all%STF,2)
do i=1, size(bc%STF,2)
call gatherv_all_cr(bc%STF(:,i),bc%dataXZ%npoin,bc%dataXZ_all%STF(:,i),bc%npoin_perproc,bc%poin_offset,&
bc%dataXZ_all%npoin,NPROC)
enddo
endif


end subroutine gather_dataXZ

Expand Down Expand Up @@ -2645,7 +2637,7 @@ subroutine fault_output_synchronize_GPU(it)
! output data
!call SCEC_write_dataT(bc%dataT)
if ( it == bc%dataT%nt ) call SCEC_write_dataT(bc%dataT)
if ( it == bc%dataT%nt .and. myrank == 0) call write_STF_GPU(bc%dataXZ_all%STF,bc%NT,bc%dataXZ_all%npoin)
if ( it == bc%NT ) call write_STF_GPU(bc%B,bc%mubulk1,bc%mubulk2,bc%STF,bc%NT,bc%nglob,myrank)
if (myrank == 0) call write_dataXZ(bc%dataXZ_all,it,ifault)
enddo

Expand Down

0 comments on commit 77abef3

Please sign in to comment.