From 77abef32f22efb4b5e894df8c2b3cb49f406a1ca Mon Sep 17 00:00:00 2001 From: elifo Date: Thu, 30 Jun 2022 21:42:30 -0700 Subject: [PATCH] write out moment rate function as stf verified in homogeneous ridgecrest model --- .../create_regions_mesh.f90 | 5 +- .../fault_generate_databases.f90 | 78 +++++++++++++++++-- src/specfem3D/fault_solver_common.f90 | 52 ++++++++----- src/specfem3D/fault_solver_dynamic.f90 | 18 ++--- 4 files changed, 114 insertions(+), 39 deletions(-) diff --git a/src/generate_databases/create_regions_mesh.f90 b/src/generate_databases/create_regions_mesh.f90 index 9db672cc3..c5d29eefa 100644 --- a/src/generate_databases/create_regions_mesh.f90 +++ b/src/generate_databases/create_regions_mesh.f90 @@ -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 @@ -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 diff --git a/src/generate_databases/fault_generate_databases.f90 b/src/generate_databases/fault_generate_databases.f90 index 70d64a3f0..bb9c4511f 100644 --- a/src/generate_databases/fault_generate_databases.f90 +++ b/src/generate_databases/fault_generate_databases.f90 @@ -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 @@ -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 @@ -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') @@ -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 @@ -555,6 +560,8 @@ subroutine save_fault_xyzcoord_ibulk(fdb) fdb%zcoordbulk2(i) = zstore_unique(K2) enddo + enddo + end subroutine save_fault_xyzcoord_ibulk @@ -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 !================================================================================= @@ -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 !==================================================================================== @@ -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 diff --git a/src/specfem3D/fault_solver_common.f90 b/src/specfem3D/fault_solver_common.f90 index 62ce9b6a2..59a1af18f 100644 --- a/src/specfem3D/fault_solver_common.f90 +++ b/src/specfem3D/fault_solver_common.f90 @@ -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: @@ -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 @@ -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 @@ -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 !------------------------------------------------------------------------ diff --git a/src/specfem3D/fault_solver_dynamic.f90 b/src/specfem3D/fault_solver_dynamic.f90 index 2a2e4c471..5d54e6a33 100644 --- a/src/specfem3D/fault_solver_dynamic.f90 +++ b/src/specfem3D/fault_solver_dynamic.f90 @@ -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') @@ -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,.. @@ -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) @@ -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 @@ -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