diff --git a/EXAMPLES/applications/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/DATA/Par_file b/EXAMPLES/applications/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/DATA/Par_file index 01a5e8972..17a57c82c 100644 --- a/EXAMPLES/applications/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/DATA/Par_file +++ b/EXAMPLES/applications/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/DATA/Par_file @@ -324,6 +324,13 @@ TRACTION_PATH = DATA/AxiSEM_tractions/3/ FKMODEL_FILE = FKmodel RECIPROCITY_AND_KH_INTEGRAL = .false. # does not work yet +#----------------------------------------------------------- +## +## Prescribed wavefield discontinuity on an interface +## +##----------------------------------------------------------- +IS_WAVEFIELD_DISCONTINUITY = .false. + #----------------------------------------------------------- # # Run modes diff --git a/EXAMPLES/applications/wavefield_discontinuity/README.md b/EXAMPLES/applications/wavefield_discontinuity/README.md new file mode 100644 index 000000000..0e4a38196 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/README.md @@ -0,0 +1,10 @@ +## Examples for wavefield injection + +**Step 1:** configure and compile SPECFEM + +**Step 2:** compile FK simulator +`cd fk_coupling; bash make_fk_injection.sh` + +**Step 3:** go inside an example directory, see a self-explanatory script `run_this_example.sh` + +**Step 4:** after running the program, plot the waveforms using `python plot_seismograms.py` diff --git a/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/compute_fk_injection_field.f90 b/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/compute_fk_injection_field.f90 new file mode 100644 index 000000000..1ea41584a --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/compute_fk_injection_field.f90 @@ -0,0 +1,254 @@ +program write_injection_field + use mpi + use fk_injection + implicit none + integer, parameter :: THREE = 3 + integer, parameter :: NGLLMID = (NGLLSQUARE + 1) / 2 + integer :: ier, ip, ib, igll, ilayer,nlines, nb, np + character(len=256) :: fn, line + real(kind=CUSTOM_REAL), dimension(:), allocatable :: xp, yp, zp + real(kind=CUSTOM_REAL), dimension(:), allocatable :: xb, yb, zb, nxb, nyb, nzb + !real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: displ, accel + !real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: traction + real(kind=CUSTOM_REAL) :: ray_p,Tg,DF_FK + + real(kind=CUSTOM_REAL) :: rho_tmp,kappa_tmp,mu_tmp,alpha_tmp,beta_tmp,xi + integer :: ii, kk, iim1, iip1, iip2, it_tmp + real(kind=CUSTOM_REAL) :: cs1,cs2,cs3,cs4,w + + real(kind=CUSTOM_REAL), parameter :: TOL_ZERO_TAKEOFF = 1.e-14 + + real(kind=CUSTOM_REAL) :: zmid, ztop + real(kind=CUSTOM_REAL) :: time_t + + call MPI_Init(ier) + if (ier /= 0 ) stop 'Error initializing MPI' + call MPI_Comm_rank(MPI_COMM_WORLD, myrank, ier) + + write(fn, "('DATABASES_MPI/proc',i6.6,'_wavefield_discontinuity_points')")& + myrank + open(77, file=trim(fn), action="read") + ier = 0 + nlines = 0 + do while (ier == 0) + read(77, '(a)', iostat=ier) line + if (ier == 0) nlines = nlines + 1 + enddo + close(77) + np = nlines + allocate(xp(np), yp(np), zp(np)) + !allocate(displ(THREE, np), accel(THREE, np)) + !displ(:,:) = 0.0 + !accel(:,:) = 0.0 + open(77, file=trim(fn), action="read") + do ip = 1, np + read(77, *) xp(ip), yp(ip), zp(ip) + enddo + close(77) + write(fn, "('DATABASES_MPI/proc',i6.6,'_wavefield_discontinuity_faces')")& + myrank + open(77, file=trim(fn), action="read") + ier = 0 + nlines = 0 + do while (ier == 0) + read(77, '(a)', iostat=ier) line + if (ier == 0) nlines = nlines + 1 + enddo + close(77) + nb = nlines / NGLLSQUARE + allocate(xb(nb*NGLLSQUARE), yb(nb*NGLLSQUARE), & + zb(nb*NGLLSQUARE), nxb(nb*NGLLSQUARE), nyb(nb*NGLLSQUARE), & + nzb(nb*NGLLSQUARE), xi1(nb*NGLLSQUARE), xim(nb*NGLLSQUARE), & + bdlambdamu(nb*NGLLSQUARE)) + !allocate(traction(THREE, NGLLSQUARE, nb)) + !traction(:,:,:) = 0.0 + open(77, file=trim(fn), action="read") + do ip = 1, nb*NGLLSQUARE + ! do igll = 1, NGLLSQUARE + read(77, *) xb(ip), yb(ip), zb(ip), & + nxb(ip), nyb(ip), nzb(ip) + ! enddo + enddo + close(77) + + + call ReadFKModelInput() + + ! converts origin point Z to reference framework depth for FK, + ! where top of lower half-space has to be at z==0 + zz0 = zz0 - Z_REF_for_FK + + ! converts to rad + phi_FK = phi_FK * PI/180.d0 ! azimuth + theta_FK = theta_FK * PI/180.d0 ! take-off + + ! ray parameter p (according to Snell's law: sin(theta1)/v1 == + ! sin(theta2)/v2) + if (type_kpsv_fk == 1) then + ! P-wave + ray_p = sin(theta_FK)/alpha_FK(nlayer) ! for vp (i.e., alpha) + else if (type_kpsv_fk == 2) then + ! SV-wave + ray_p = sin(theta_FK)/beta_FK(nlayer) ! for vs (i.e., beta) + endif + + ! note: vertical incident (theta==0 -> p==0) is not handled. + ! here, it limits ray parameter p to a very small value to handle + ! the calculations + if (abs(ray_p) < TOL_ZERO_TAKEOFF) ray_p = sign(TOL_ZERO_TAKEOFF,ray_p) + + ! maximum period + Tg = 1.d0 / ff0 + + call find_size_of_working_arrays(deltat, freq_sampling_fk, tmax_fk, & + NF_FOR_STORING, & + NF_FOR_FFT, NPOW_FOR_INTERP, NP_RESAMP, & + DF_FK) + + ! user output + if (myrank == 0) then + write(IMAIN,*) ' computed FK parameters:' + write(IMAIN,*) ' frequency sampling rate = ', freq_sampling_fk,"(Hz)" + write(IMAIN,*) ' number of frequencies to store = ', NF_FOR_STORING + write(IMAIN,*) ' number of frequencies for FFT = ', NF_FOR_FFT + write(IMAIN,*) ' power of 2 for FFT = ', NPOW_FOR_INTERP + write(IMAIN,*) + write(IMAIN,*) ' simulation time step = ', deltat,"(s)" + write(IMAIN,*) ' total simulation length = ', NSTEP*deltat,"(s)" + write(IMAIN,*) + write(IMAIN,*) ' FK time resampling rate = ', NP_RESAMP + write(IMAIN,*) ' new time step for F-K = ', NP_RESAMP * deltat,"(s)" + write(IMAIN,*) ' new time window length = ', tmax_fk,"(s)" + write(IMAIN,*) + write(IMAIN,*) ' frequency step for F-K = ', DF_FK,"(Hz)" + call flush_IMAIN() + endif + + ! safety check with number of simulation time steps + if (NSTEP/NP_RESAMP > NF_FOR_STORING + NP_RESAMP) then + if (myrank == 0) then + print *,'Error: FK time window length ',tmax_fk,' and NF_for_storing ',NF_FOR_STORING + print *,' are too small for chosen simulation length with NSTEP = ',NSTEP + print * + print *,' you could use a smaller NSTEP <= ',NF_FOR_STORING*NP_RESAMP + print *,' or' + print *,' increase FK window length larger than ',(NSTEP/NP_RESAMP - NP_RESAMP) * NP_RESAMP * deltat + print *,' to have a NF for storing larger than ',(NSTEP/NP_RESAMP - NP_RESAMP) + endif + stop 'Invalid FK setting' + endif + + ! safety check + if (NP_RESAMP == 0) then + if (myrank == 0) then + print *,'Error: FK resampling rate ',NP_RESAMP,' is invalid for frequency sampling rate ',freq_sampling_fk + print *,' and the chosen simulation DT = ',deltat + print * + print *,' you could use a higher frequency sampling rate>',1./(deltat) + print *,' (or increase the time stepping size DT if possible)' + endif + stop 'Invalid FK setting' + endif + + ! limits resampling sizes + if (NP_RESAMP > 10000) then + if (myrank == 0) then + print *,'Error: FK resampling rate ',NP_RESAMP,' is too high for frequency sampling rate ',freq_sampling_fk + print *,' and the chosen simulation DT = ',deltat + print * + print *,' you could use a higher frequency sampling rate>',1./(10000*deltat) + print *,' (or increase the time stepping size DT if possible)' + endif + stop 'Invalid FK setting' + endif + + allocate(displ(THREE, np, -NP_RESAMP:NF_FOR_STORING+NP_RESAMP), & + accel(THREE, np, -NP_RESAMP:NF_FOR_STORING+NP_RESAMP), & + traction(THREE, NGLLSQUARE*nb, -NP_RESAMP:NF_FOR_STORING+NP_RESAMP)) + + do ib = 1, nb + zmid = zb((ib-1)*NGLLSQUARE+NGLLMID) + ztop = 0.0 + do ilayer = 1, nlayer-1 + ztop = ztop - h_FK(ilayer) + if (ztop < zmid) exit + enddo + if (ztop > zmid) ilayer = nlayer + rho_tmp = rho_FK(ilayer) + alpha_tmp = alpha_FK(ilayer) + beta_tmp = beta_FK(ilayer) + kappa_tmp = rho_tmp*(alpha_tmp*alpha_tmp-4.0/3.0*beta_tmp*beta_tmp) + mu_tmp = rho_tmp*beta_tmp*beta_tmp + xi = mu_tmp/(kappa_tmp + 4.0/3.0 * mu_tmp) + xi1((ib-1)*NGLLSQUARE+1:ib*NGLLSQUARE) = 1.0 - 2.0 * xi + xim((ib-1)*NGLLSQUARE+1:ib*NGLLSQUARE) = (1.0 - xi) * mu_tmp + bdlambdamu((ib-1)*NGLLSQUARE+1:ib*NGLLSQUARE) = & + (3.0 * kappa_tmp - 2.0 * mu_tmp) / (6.0 * kappa_tmp + 2.0 * mu_tmp) ! Poisson's ratio 3K-2G/[2(3K+G)] + enddo + + allocate(xx(np), yy(np), zz(np)) + xx(:) = xp(:); yy(:) = yp(:); zz(:) = zp(:) - Z_REF_for_FK + call FK(alpha_FK, beta_FK, mu_FK, h_FK, nlayer, & + Tg, ray_p, phi_FK, xx0, yy0, zz0, & + tt0, deltat, nstep, np, & + type_kpsv_fk, NF_FOR_STORING, NPOW_FOR_FFT, NP_RESAMP, DF_FK, & + .false.) + deallocate(xx, yy, zz) + allocate(xx(nb*NGLLSQUARE), yy(nb*NGLLSQUARE), zz(nb*NGLLSQUARE), & + nmx(nb*NGLLSQUARE), nmy(nb*NGLLSQUARE), nmz(nb*NGLLSQUARE)) + xx(:) = xb(:); yy(:) = yb(:); zz(:) = zb(:) - Z_REF_for_FK + nmx(:) = nxb(:); nmy(:) = nyb(:); nmz(:) = nzb(:) + call FK(alpha_FK, beta_FK, mu_FK, h_FK, nlayer, & + Tg, ray_p, phi_FK, xx0, yy0, zz0, & + tt0, deltat, nstep, nb*NGLLSQUARE, & + type_kpsv_fk, NF_FOR_STORING, NPOW_FOR_FFT, NP_RESAMP, DF_FK, & + .true.) + deallocate(xx, yy, zz) + write(fn, "('DATABASES_MPI/proc',i6.6,'_wavefield_discontinuity.bin')")& + myrank + open(88, file=trim(fn), form="unformatted", action="write") + if (myrank == 0) open(99, file='injection_displ', form='formatted', action='write') + ! data + do it_tmp = 1,NSTEP + ! FK coupling + !! find indices + ! example: + ! np_resamp = 1 and it = 1,2,3,4,5,6, .. + ! --> ii = 1,2,3,4,5,6,.. + ! np_resamp = 2 and it = 1,2,3,4,5,6, .. + ! --> ii = 1,1,2,2,3,3,.. + ii = floor( real(it_tmp + NP_RESAMP - 1) / real( NP_RESAMP)) + ! example: + ! kk = 1,2,1,2,1,2,,.. + kk = it_tmp - (ii-1) * NP_RESAMP + ! example: + ! w = 0,1/2,0,1/2,.. + w = dble(kk-1) / dble(NP_RESAMP) + + ! Cubic spline values + cs4 = w*w*w/6.d0 + cs1 = 1.d0/6.d0 + w*(w-1.d0)/2.d0 - cs4 + cs3 = w + cs1 - 2.d0*cs4 + cs2 = 1.d0 - cs1 - cs3 - cs4 + + ! interpolation indices + iim1 = ii-1 ! 0,.. + iip1 = ii+1 ! 2,.. + iip2 = ii+2 ! 3,.. + + write(88) cs1*displ(:,:,iim1)+cs2*displ(:,:,ii)+cs3*displ(:,:,iip1)& + +cs4*displ(:,:,iip2) + write(88) cs1*accel(:,:,iim1)+cs2*accel(:,:,ii)+cs3*accel(:,:,iip1)& + +cs4*accel(:,:,iip2) + write(88) cs1*traction(:,:,iim1)+cs2*traction(:,:,ii)& + +cs3*traction(:,:,iip1)+cs4*traction(:,:,iip2) + ! time + time_t = (it_tmp-1) * deltat - tt0 + if (myrank == 0) write(99, *) time_t, & + cs1*displ(3,1,iim1)+cs2*displ(3,1,ii)+cs3*displ(3,1,iip1)& + +cs4*displ(3,1,iip2) + enddo + close(88) + close(99) + call MPI_Finalize(ier) +end program write_injection_field diff --git a/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/compute_fk_receiver.f90 b/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/compute_fk_receiver.f90 new file mode 100644 index 000000000..8ee15f5e6 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/compute_fk_receiver.f90 @@ -0,0 +1,174 @@ +program compute_fk_receiver + use fk_injection + implicit none + integer, parameter :: THREE = 3 + character(len=256) :: fn, line + character(len=3) :: char_comp = 'XYZ' + integer :: ier, ip, np, nlines, ic + real(kind=CUSTOM_REAL), dimension(:), allocatable :: xp, yp, zp + real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: seismogram + character(len=256), dimension(:), allocatable :: nt_name, sta_name + real(kind=CUSTOM_REAL) :: ray_p,Tg,DF_FK + + real(kind=CUSTOM_REAL) :: rho_tmp,kappa_tmp,mu_tmp,alpha_tmp,beta_tmp,xi + integer :: ii, kk, iim1, iip1, iip2, it_tmp + real(kind=CUSTOM_REAL) :: cs1,cs2,cs3,cs4,w + + real(kind=CUSTOM_REAL), parameter :: TOL_ZERO_TAKEOFF = 1.e-14 + + real(kind=CUSTOM_REAL) :: zmid, ztop + real(kind=CUSTOM_REAL) :: time_t + + fn = "DATA/STATIONS" + open(77, file=trim(fn), action="read") + ier = 0 + nlines = 0 + do while (ier == 0) + read(77, '(a)', iostat=ier) line + if (ier == 0) nlines = nlines + 1 + enddo + close(77) + np = nlines + allocate(xp(np), yp(np), zp(np), nt_name(np), sta_name(np)) + open(77, file=trim(fn), action="read") + do ip = 1, np + read(77, *) sta_name(ip), nt_name(ip), yp(ip), xp(ip), zp(ip) + enddo + close(77) + + call ReadFKModelInput() + + ! converts origin point Z to reference framework depth for FK, + ! where top of lower half-space has to be at z==0 + zz0 = zz0 - Z_REF_for_FK + + ! converts to rad + phi_FK = phi_FK * PI/180.d0 ! azimuth + theta_FK = theta_FK * PI/180.d0 ! take-off + + ! ray parameter p (according to Snell's law: sin(theta1)/v1 == + ! sin(theta2)/v2) + if (type_kpsv_fk == 1) then + ! P-wave + ray_p = sin(theta_FK)/alpha_FK(nlayer) ! for vp (i.e., alpha) + else if (type_kpsv_fk == 2) then + ! SV-wave + ray_p = sin(theta_FK)/beta_FK(nlayer) ! for vs (i.e., beta) + endif + + ! note: vertical incident (theta==0 -> p==0) is not handled. + ! here, it limits ray parameter p to a very small value to handle + ! the calculations + if (abs(ray_p) < TOL_ZERO_TAKEOFF) ray_p = sign(TOL_ZERO_TAKEOFF,ray_p) + + ! maximum period + Tg = 1.d0 / ff0 + + call find_size_of_working_arrays(deltat, freq_sampling_fk, tmax_fk, & + NF_FOR_STORING, & + NF_FOR_FFT, NPOW_FOR_INTERP, NP_RESAMP, & + DF_FK) + + ! user output + if (myrank == 0) then + write(IMAIN,*) ' computed FK parameters:' + write(IMAIN,*) ' frequency sampling rate = ', freq_sampling_fk,"(Hz)" + write(IMAIN,*) ' number of frequencies to store = ', NF_FOR_STORING + write(IMAIN,*) ' number of frequencies for FFT = ', NF_FOR_FFT + write(IMAIN,*) ' power of 2 for FFT = ', NPOW_FOR_INTERP + write(IMAIN,*) + write(IMAIN,*) ' simulation time step = ', deltat,"(s)" + write(IMAIN,*) ' total simulation length = ', NSTEP*deltat,"(s)" + write(IMAIN,*) + write(IMAIN,*) ' FK time resampling rate = ', NP_RESAMP + write(IMAIN,*) ' new time step for F-K = ', NP_RESAMP * deltat,"(s)" + write(IMAIN,*) ' new time window length = ', tmax_fk,"(s)" + write(IMAIN,*) + write(IMAIN,*) ' frequency step for F-K = ', DF_FK,"(Hz)" + call flush_IMAIN() + endif + + ! safety check with number of simulation time steps + if (NSTEP/NP_RESAMP > NF_FOR_STORING + NP_RESAMP) then + if (myrank == 0) then + print *,'Error: FK time window length ',tmax_fk,' and NF_for_storing ',NF_FOR_STORING + print *,' are too small for chosen simulation length with NSTEP = ',NSTEP + print * + print *,' you could use a smaller NSTEP <= ',NF_FOR_STORING*NP_RESAMP + print *,' or' + print *,' increase FK window length larger than ',(NSTEP/NP_RESAMP -NP_RESAMP) * NP_RESAMP * deltat + print *,' to have a NF for storing larger than ',(NSTEP/NP_RESAMP -NP_RESAMP) + endif + stop 'Invalid FK setting' + endif + + ! limits resampling sizes + if (NP_RESAMP > 10000) then + if (myrank == 0) then + print *,'Error: FK resampling rate ',NP_RESAMP,' is too high for frequencysampling rate ',freq_sampling_fk + print *,' and the chosen simulation DT = ',deltat + print * + print *,' you could use a higher frequency sampling rate>',1./(10000*deltat) + print *,' (or increase the time stepping size DT if possible)' + endif + stop 'Invalid FK setting' + endif + + allocate(displ(THREE, np, -NP_RESAMP:NF_FOR_STORING+NP_RESAMP), & + accel(THREE, np, -NP_RESAMP:NF_FOR_STORING+NP_RESAMP)) + + allocate(xx(np), yy(np), zz(np)) + xx(:) = xp(:); yy(:) = yp(:); zz(:) = zp(:) - Z_REF_for_FK + call FK(alpha_FK, beta_FK, mu_FK, h_FK, nlayer, & + Tg, ray_p, phi_FK, xx0, yy0, zz0, & + tt0, deltat, nstep, np, & + type_kpsv_fk, NF_FOR_STORING, NPOW_FOR_FFT, NP_RESAMP, DF_FK, & + .false.) + deallocate(xx, yy, zz) + + allocate(seismogram(THREE, np, NSTEP)) + do it_tmp = 1,NSTEP + ! FK coupling + !! find indices + ! example: + ! np_resamp = 1 and it = 1,2,3,4,5,6, .. + ! --> ii = 1,2,3,4,5,6,.. + ! np_resamp = 2 and it = 1,2,3,4,5,6, .. + ! --> ii = 1,1,2,2,3,3,.. + ii = floor( real(it_tmp + NP_RESAMP - 1) / real( NP_RESAMP)) + ! example: + ! kk = 1,2,1,2,1,2,,.. + kk = it_tmp - (ii-1) * NP_RESAMP + ! example: + ! w = 0,1/2,0,1/2,.. + w = dble(kk-1) / dble(NP_RESAMP) + + ! Cubic spline values + cs4 = w*w*w/6.d0 + cs1 = 1.d0/6.d0 + w*(w-1.d0)/2.d0 - cs4 + cs3 = w + cs1 - 2.d0*cs4 + cs2 = 1.d0 - cs1 - cs3 - cs4 + + ! interpolation indices + iim1 = ii-1 ! 0,.. + iip1 = ii+1 ! 2,.. + iip2 = ii+2 ! 3,.. + + seismogram(:,:,it_tmp) = cs1*displ(:,:,iim1)+cs2*displ(:,:,ii)+& + cs3*displ(:,:,iip1)+cs4*displ(:,:,iip2) + enddo + + do ip = 1, np + do ic = 1, THREE + fn = "OUTPUT_FILES/"//trim(nt_name(ip))//'.'//trim(sta_name(ip)) & + //'.'//char_comp(ic:ic)//".fkd" + open(77, file=trim(fn), action='write') + do it_tmp = 1, NSTEP + time_t = (it_tmp-1) * deltat - tt0 + write(77, *) time_t, seismogram(ic, ip, it_tmp) + enddo + close(77) + enddo + enddo + deallocate(xp, yp, zp, nt_name, sta_name, displ, accel, seismogram) +end program compute_fk_receiver diff --git a/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/coupling_fk.f90 b/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/coupling_fk.f90 new file mode 100644 index 000000000..9720150d7 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/coupling_fk.f90 @@ -0,0 +1,1529 @@ + + +module fk_injection + integer, parameter :: CUSTOM_REAL = 4, IMAIN = 6, NGLLSQUARE = 25 + double precision, parameter :: PI = 3.141592653589793d0, & + TINYVAL = 1.d-9 + integer :: myrank + integer :: NSTEP + real(kind=CUSTOM_REAL) :: deltat + character(len=100) :: FKMODEL_FILE = 'fk_model' + ! FK elastic + integer :: nlayer + integer :: NF_FOR_STORING, NF_FOR_FFT, NPOW_FOR_FFT, NP_RESAMP, NPOW_FOR_INTERP + integer :: NPTS_STORED, NPTS_INTERP + + !real(kind=CUSTOM_REAL),dimension(:,:),allocatable :: & + ! vxbd,vybd,vzbd,txxbd,txybd,txzbd,tyybd,tyzbd,tzzbd + real(kind=CUSTOM_REAL) :: Z_REF_for_FK + + ! source + integer :: type_kpsv_fk = 0 ! incident wave type: 1 == P-wave, 2 == SV-wave + real(kind=CUSTOM_REAL) :: xx0,yy0,zz0,ff0,tt0,tmax_fk,freq_sampling_fk,amplitude_fk + real(kind=CUSTOM_REAL) :: phi_FK,theta_FK + + ! model + real(kind=CUSTOM_REAL),dimension(:),allocatable :: alpha_FK,beta_FK,rho_FK,mu_FK,h_FK + !complex(kind=8), dimension(:,:), allocatable :: VX_f, VY_f, VZ_f, TX_f, TY_f, TZ_f + real(kind=CUSTOM_REAL),dimension(:,:,:),allocatable :: displ, accel, traction + + !complex(kind=8), dimension(:), allocatable :: WKS_CMPLX_FOR_FFT + !real(kind=CUSTOM_REAL), dimension(:), allocatable :: WKS_REAL_FOR_FFT + + real(kind=CUSTOM_REAL),dimension(:),allocatable :: xx,yy,zz,xi1,xim,bdlambdamu + + ! normal + real(kind=CUSTOM_REAL),dimension(:),allocatable :: nmx,nmy,nmz +end module fk_injection + + subroutine FK(al, be, mu, H, nlayer, & + Tg, ray_p, phi, x0, y0, z0, & + t0, dt, npts, npt, & + kpsv, NF_FOR_STORING, NPOW_FOR_FFT, NP_RESAMP, DF_FK, & + compute_traction) + + use fk_injection, only: myrank,CUSTOM_REAL,IMAIN,PI,TINYVAL,& + NGLLSQUARE, & + displ, accel, traction, & + xx, yy, zz, xi1, xim, bdlambdamu, & + nmx, nmy, nmz, NPTS_STORED, NPTS_INTERP, & + amplitude_fk + + implicit none + + integer, parameter :: CUSTOM_CMPLX = 8 + real(kind=CUSTOM_REAL), parameter :: zign_neg = -1.0 + + ! input and output + integer, intent(in) :: nlayer, npt, npts, kpsv + integer :: NF_FOR_STORING, NPOW_FOR_FFT, NP_RESAMP + + ! model + real(kind=CUSTOM_REAL), dimension(nlayer), intent(in) :: al(nlayer),be(nlayer),mu(nlayer),H(nlayer) + + ! source + real(kind=CUSTOM_REAL), intent(in) :: dt, ray_p, phi, x0, y0, z0, Tg, t0, DF_FK + + ! local parameters + real(kind=CUSTOM_REAL), dimension(:), allocatable :: fvec, dtmp + real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: field + real(kind=CUSTOM_REAL), dimension(:), allocatable :: tmp_t1, tmp_t2, tmp_t3, tmp_it1 + + complex(kind=CUSTOM_CMPLX), dimension(:,:), allocatable :: coeff, field_f + complex(kind=CUSTOM_CMPLX), dimension(:), allocatable :: tmp_f1, tmp_f2, tmp_f3 + complex(kind=CUSTOM_CMPLX) :: C_3,stf_coeff,a,b,c,d,delta_mat + complex(kind=CUSTOM_CMPLX) :: dx_f,dz_f,txz_f,tzz_f + complex(kind=CUSTOM_CMPLX) :: N_mat(4,4) + + real(kind=CUSTOM_REAL) :: sigma_rr,sigma_rt,sigma_rz,sigma_tt,sigma_tz,sigma_zz + real(kind=CUSTOM_REAL) :: Txx_tmp, Txy_tmp, Txz_tmp, Tyy_tmp, Tyz_tmp, Tzz_tmp + real(kind=CUSTOM_REAL) :: dt_fk,df,om,Tdelay,eta_p,eta_s,f_Nyquist,C_1 + + integer :: npow,npts2,nf,nf2,nn,ii,ipt,i,j,lpts + integer :: npoints2 + integer :: ier,iface,igll,ispec + + ! spline work array + double precision, dimension(:), allocatable :: tmp_c + +!! DK DK here is the hardwired maximum size of the array +!! DK DK Aug 2016: if this routine is called many times (for different mesh points at which the SEM is coupled with FK) +!! DK DK Aug 2016: this should be moved to the calling program and precomputed once and for all + real(kind=CUSTOM_REAL) :: mpow(30) + + ! taper + ! idea: tapers the onset of the injection traces to diminish numerical noise. + ! the inverse FFT can lead to non-zero onsets for coarse frequency sampling. + logical, parameter :: USE_TAPERED_BEGINNING = .false. + integer, parameter :: taper_nlength = 20 ! tapers first 20 steps + real(kind=CUSTOM_REAL) :: taper + + ! fixed parameters + !integer, parameter :: nvar = 5 + logical :: compute_traction + integer :: nvar + + if (compute_traction) then + nvar = 3 + else + nvar = 4 + endif + + ! initializations + !! new way to do time domain resampling + df = DF_FK + nf2 = NF_FOR_STORING+1 ! number of positive frequency sample points + nf = 2*NF_FOR_STORING ! number of total frequencies after symmetrisation + npts2 = nf ! number of samples in time serie + + !! VM VM recompute new values for new way to do + npow = ceiling(log(npts2*1.0)/log(2.0)) + npts2 = 2**npow + NPOW_FOR_FFT = npow + + !! number of points for resampled vector + npoints2 = NP_RESAMP*(npts2-1)+1 + + dt_fk = 1.0_CUSTOM_REAL/(df*(npts2-1)) + f_Nyquist = 1.0_CUSTOM_REAL/(2.0_CUSTOM_REAL * dt) ! Nyquist frequency of specfem time serie + + ! user output + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) ' Entering the FK synthetics program:' + write(IMAIN,*) ' Number of samples stored for FK solution = ', NF_FOR_STORING + write(IMAIN,*) ' Number of points used for FFT = ', npts2 + write(IMAIN,*) ' Total time length used for FK = ', t0+(npts2-1)*dt_fk,'(s)' + write(IMAIN,*) + write(IMAIN,*) ' simulation Nyquist frequency = ', f_Nyquist,'(Hz)' + write(IMAIN,*) ' FK time step = ', dt_fk + write(IMAIN,*) ' FK frequency step = ', df + write(IMAIN,*) ' power of 2 for FFT = ', npow + write(IMAIN,*) + call flush_IMAIN() + endif + + !! check if dt_fk is compatible with dt_specfem + !! + !! + !! + + allocate(fvec(nf2),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating array 2218') + if (ier /= 0) stop 'error while allocating' + fvec(:) = 0.0_CUSTOM_REAL + do ii = 1, nf2 + fvec(ii) = (ii-1)*df + enddo + + allocate(coeff(2,nf2),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating array 2219') + if (ier /= 0) stop 'error while allocating' + coeff(:,:) = (0.d0,0.d0) + + allocate(field_f(nf,nvar),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating array 2220') + if (ier /= 0) stop 'error while allocating' + field_f(:,:) = (0.d0,0.d0) + + allocate(field(npts2,nvar),dtmp(npts),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating array 2221') + if (ier /= 0) stop 'error while allocating' + field(:,:) = 0._CUSTOM_REAL; dtmp(:) = 0._CUSTOM_REAL + + !! allocate debug vectors + allocate(tmp_f1(npts2), tmp_f2(npts2), tmp_f3(npts2),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating array 2222') + if (ier /= 0) stop 'error while allocating' + tmp_f1(:) = (0.d0,0.d0) + tmp_f2(:) = (0.d0,0.d0) + tmp_f3(:) = (0.d0,0.d0) + + if (ier /= 0) stop 'error while allocating' + allocate(tmp_t1(npts2), tmp_t2(npts2), tmp_t3(npts2),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating array 2223') + if (ier /= 0) stop 'error while allocating' + tmp_t1(:) = 0.0_CUSTOM_REAL + tmp_t2(:) = 0.0_CUSTOM_REAL + tmp_t3(:) = 0.0_CUSTOM_REAL + + allocate(tmp_it1(npoints2),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating array 2224') + if (ier /= 0) stop 'error while allocating' + tmp_it1(:) = 0.0_CUSTOM_REAL + + ! temporary work array for splines + allocate(tmp_c(npts2),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating array 2225') + tmp_c(:) = 0.d0 + + NPTS_STORED = npts2 + NPTS_INTERP = npoints2 + + nn = int(-t0/dt) ! what if this is not an integer number? + +!! DK DK Aug 2016: if this routine is called many times (for different mesh points at which the SEM is coupled with FK) +!! DK DK Aug 2016: this should be moved to the calling program and precomputed once and for all + do i = 1,npow + mpow(i) = 2**(npow-i) + enddo + + if (myrank == 0) then + write(IMAIN,*) ' starting from ',nn,' points before time 0' + write(IMAIN,*) + call flush_IMAIN() + endif + + if (kpsv == 1) then + ! P-wave + ! for C_3 = i sin(inc) (u=[sin(inc), cos(inc)]) + C_3 = amplitude_fk * cmplx(0,1.) * ray_p * al(nlayer) ! amp. of incoming P in the bot. layer + eta_p = sqrt(1.0/al(nlayer)**2 - ray_p**2) ! vertical slowness for lower layer + + if (myrank == 0) write(IMAIN,*) ' Incoming P : C_3, ray_p, eta = ', C_3, ray_p, eta_p + + N_mat(:,:) = (0.d0,0.d0) + + ! find out the wave coefficients in the bottom layer for all freqs ------------------------------- + do ii = 1, nf2 + om = 2.0 * PI * fvec(ii) + ! propagation matrix + call compute_N_Rayleigh(al,be,mu,H,nlayer,om,ray_p,sum(H(1:nlayer-1)),N_mat) !total-thickness=sum(H) + + a = N_mat(3,2); b = N_mat(3,4); c = N_mat(4,2); d = N_mat(4,4) + delta_mat = a*d - b*c + if (abs(delta_mat) > TINYVAL) then + coeff(1,ii) = -(d*N_mat(3,3) - b*N_mat(4,3)) / delta_mat * C_3 + coeff(2,ii) = -(-c*N_mat(3,3) + a*N_mat(4,3)) / delta_mat * C_3 + else + coeff(1,ii) = (0.d0,0.d0) + coeff(2,ii) = (0.d0,0.d0) + endif + + !debug + !if (ii == 1 .and. myrank == 0) then + ! print *,'debug: Rayleigh coeff ',coeff(1,ii),coeff(2,ii),delta_mat + ! print *,'N_mat' + ! print *,N_mat + !endif + enddo + + ! loop over all data points ------------------------------------------------- + ! instead of looping over points like: + do ipt = 1, npt ! maybe this can be run faster by shifting t for diff. x of fixed z + ! we loop over the boundary arrays to get ispec & acoustic/elastic domain flag: + !do iface = 1,num_abs_boundary_faces + ! ispec = abs_boundary_ispec(iface) + + ! ! GLL points on boundary face + ! do igll = 1,NGLLSQUARE + ! ! point index using table lookup + ! ipt = ipt_table(igll,iface) + + ! initializes + field_f(:,:) = (0.d0,0.d0) + + ! time delay with respect to top of lower half-space (set to be at z==0) + Tdelay = ray_p * (xx(ipt)-x0) * cos(phi) + ray_p * (yy(ipt)-y0) * sin(phi) + eta_p * (0.0-z0) + + do ii = 1, nf2 + om = 2.0 * PI * fvec(ii) !! pulsation + + stf_coeff = exp(-(om * Tg/2)**2) !! apodization window + stf_coeff = stf_coeff * exp(cmplx(0,-1)*om*Tdelay) + + !! zz(ipt) is the height of point with respect to the lower layer + call compute_N_Rayleigh(al,be,mu,H,nlayer,om,ray_p,zz(ipt),N_mat) + + + dx_f = N_mat(1,2)*coeff(1,ii) + N_mat(1,4)*coeff(2,ii) + N_mat(1,3)*C_3 ! y_1 + dz_f = N_mat(2,2)*coeff(1,ii) + N_mat(2,4)*coeff(2,ii) + N_mat(2,3)*C_3 ! y_3 + if (.not. compute_traction) then + ! for the Stacey boundary contribution, we need velocity = (i om) displacement (in frequency domain) + field_f(ii,1) = stf_coeff * dx_f * cmplx(0,-1) ! u_x + field_f(ii,2) = stf_coeff * dz_f ! u_z + + field_f(ii,3) = stf_coeff * dx_f * cmplx(0,-1) * cmplx(-om*om,0) ! (-om^2)u_x + field_f(ii,4) = stf_coeff * dz_f * cmplx(-om*om,0) ! (-om^2)u_z + ! acoustic boundary point + ! note: instead of velocity as in the elastic case, in acoustic domains we would need potentials. + ! the velocity potential would be defined as: v = 1/rho grad(potential_dot) + ! thus, we would require to either change the FK formulations or to integrate velocity. + ! this will be left as a todo for future... + ! + ! as a reminder, displacement in frequency domains could be obtained by: + !if (ispec_is_acoustic(ispec)) then + ! ! displacement: u_x = -i y_1 + ! ! u_z = y_3 from Tong. (A13) + ! field_f(ii,1) = stf_coeff * dx_f * cmplx(0,-1) ! u_x = - i y_1 + ! field_f(ii,2) = stf_coeff * dz_f ! u_z = y_3 + !endif + + ! stress + else + txz_f = N_mat(3,2)*coeff(1,ii) + N_mat(3,4)*coeff(2,ii) + N_mat(3,3)*C_3 ! tilde{y}_4 + tzz_f = N_mat(4,2)*coeff(1,ii) + N_mat(4,4)*coeff(2,ii) + N_mat(4,3)*C_3 ! tilde{y}_6 + + field_f(ii,1) = stf_coeff * om * ray_p * (xi1(ipt)*tzz_f - 4.0*xim(ipt)*dx_f) ! T_xx + field_f(ii,2) = stf_coeff * om * ray_p * txz_f * cmplx(0,-1) ! T_xz + field_f(ii,3) = stf_coeff * om * ray_p * tzz_f ! T_zz + endif + + !debug + !if (ipt==1000 .and. ii == 10 .and. myrank == 0) print *,'debug: coeff',coeff(1,ii),coeff(2,ii), & + ! 'dx_f',dx_f,tzz_f,'xi',xi1(ipt),xim(ipt),'field',field_f(ii,3:5) + enddo + + ! pad negative f, and convert to time series + do ii = 2, nf2-1 + field_f(nf+2-ii,:) = conjg(field_f(ii,:)) + enddo + + !! inverse fast fourier transform + field(:,:) = 0.0 + do j = 1, nvar + ! inverse FFT + call FFTinv(npow,field_f(:,j),zign_neg,dt,field(:,j),mpow) + + ! wrap around to start from t0: here one has to be careful if t0/dt is not + ! exactly an integer, assume nn > 0 + if (nn > 0) then + dtmp(1:nn) = field(npts2-nn+1:npts2,j) + field(nn+1:npts2,j) = field(1:npts2-nn,j) + field(1:nn,j) = dtmp(1:nn) + else if (nn < 0) then + dtmp(1:nn) = field(1:nn,j) + field(1:npts-nn,j) = field(nn+1:npts,j) + field(npts-nn+1:npts,j) = dtmp(1:nn) + endif + enddo + + if (.not. compute_traction) then + !! store undersampled version of velocity FK solution + tmp_t1(:) = field(:,1) * cos(phi) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + displ(1,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(:) = field(:,1) * sin(phi) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + displ(2,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(:) = field(:,2) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + displ(3,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(:) = field(:,3) * cos(phi) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + accel(1,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(:) = field(:,3) * sin(phi) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + accel(2,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(:) = field(:,4) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + accel(3,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + else + !! compute traction + do lpts = 1, NF_FOR_STORING + sigma_rr = field(lpts,1) + sigma_rt = 0.0 + sigma_rz = field(lpts,2) + sigma_zz = field(lpts,3) + sigma_tt = bdlambdamu(ipt)*(sigma_rr+sigma_zz) + sigma_tz = 0.0 + + Txx_tmp = sigma_rr * cos(phi) * cos(phi) + sigma_tt * sin(phi) * sin(phi) + Txy_tmp = cos(phi) * sin(phi) * (sigma_rr - sigma_tt) + Txz_tmp = sigma_rz * cos(phi) + Tyy_tmp = sigma_rr * sin(phi) * sin(phi) + sigma_tt * cos(phi) * cos(phi) + Tyz_tmp = sigma_rz * sin(phi) + Tzz_tmp = sigma_zz + + !! store directly the traction + traction(1,ipt,lpts) = Txx_tmp*nmx(ipt) + Txy_tmp*nmy(ipt) + Txz_tmp*nmz(ipt) + traction(2,ipt,lpts) = Txy_tmp*nmx(ipt) + Tyy_tmp*nmy(ipt) + Tyz_tmp*nmz(ipt) + traction(3,ipt,lpts) = Txz_tmp*nmx(ipt) + Tyz_tmp*nmy(ipt) + Tzz_tmp*nmz(ipt) + enddo + + !! store undersamped version of tractions FK solution + tmp_t1(1:NF_FOR_STORING) = traction(1,ipt,1:NF_FOR_STORING) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + traction(1,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(1:NF_FOR_STORING) = traction(2,ipt,1:NF_FOR_STORING) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + traction(2,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(1:NF_FOR_STORING) = traction(3,ipt,1:NF_FOR_STORING) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + traction(3,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + endif + + ! user output + if (myrank == 0 .and. npt > 1000) then + if (mod(ipt,(npt/10)) == 0) then + write(IMAIN,*) ' done ',ipt/(npt/10)*10,'% points out of ',npt + call flush_IMAIN() + endif + endif + !enddo + enddo + + else if (kpsv == 2) then + ! SV-wave + ! for C_2 = sin(inc) (u=[cos(inc), sin(inc)]) + C_1 = amplitude_fk * ray_p * be(nlayer) ! amp. of incoming S in the bot. layer + eta_s = sqrt(1.0/be(nlayer)**2 - ray_p**2) ! vertical slowness for lower layer + + if (myrank == 0) write(IMAIN,*) ' Incoming S : C_1, ray_p, eta = ', C_1, ray_p, eta_s + + N_mat(:,:) = (0.d0,0.d0) + + ! find out the wave coefficients in the bottom layer for all freqs + do ii = 1, nf2 + om = 2.0 * PI * fvec(ii) + ! propagation matrix + call compute_N_Rayleigh(al,be,mu,H,nlayer,om,ray_p,sum(H(1:nlayer-1)),N_mat) !total-thickness=sum(h) + + a = N_mat(3,2); b = N_mat(3,4); c = N_mat(4,2); d = N_mat(4,4) + delta_mat = a*d - b*c + if (abs(delta_mat) > TINYVAL) then + coeff(1,ii) = -(d*N_mat(3,1) - b*N_mat(4,1)) / delta_mat * C_1 + coeff(2,ii) = -(-c*N_mat(3,1) + a*N_mat(4,1)) / delta_mat * C_1 + else + coeff(1,ii) = (0.d0,0.d0) + coeff(2,ii) = (0.d0,0.d0) + endif + enddo + + ! loop over all data points + ! instead of looping over points like: + do ipt = 1, npt ! maybe this can be run faster by shifting t for diff. x of fixed z + ! we loop over the boundary arrays to get ispec & acoustic/elastic domain flag: + !do iface = 1,num_abs_boundary_faces + ! ispec = abs_boundary_ispec(iface) + + ! ! GLL points on boundary face + ! do igll = 1,NGLLSQUARE + ! point index using table lookup + !ipt = ipt_table(igll,iface) + + ! initializes + field_f(:,:) = (0.d0,0.d0) + + ! time delay with respect to top of lower half-space (set to be at z==0) + Tdelay = ray_p * (xx(ipt)-x0) * cos(phi) + ray_p * (yy(ipt)-y0) * sin(phi) + eta_s * (0.0-z0) + + do ii = 1, nf2 + om = 2.0 * PI * fvec(ii) + stf_coeff = exp(-(om * Tg/2)**2) * exp(cmplx(0,-1)*om*Tdelay) + + ! z is the height of position with respect to the lowest layer interface. + call compute_N_Rayleigh(al,be,mu,H,nlayer,om,ray_p,zz(ipt),N_mat) + + + dx_f = N_mat(1,2)*coeff(1,ii) + N_mat(1,4)*coeff(2,ii) + N_mat(1,1)*C_1 ! y_1 + dz_f = N_mat(2,2)*coeff(1,ii) + N_mat(2,4)*coeff(2,ii) + N_mat(2,1)*C_1 ! y_3 + if (compute_traction) then + ! for the Stacey boundary contribution, we need velocity = (i om) displacement (in frequency domain) + field_f(ii,1) = stf_coeff * dx_f * cmplx(0,-1) ! u_x(1.20) + field_f(ii,2) = stf_coeff * dz_f ! u_z + field_f(ii,3) = stf_coeff * dx_f * cmplx(0,-1) * cmplx(-om*om,0) ! (-om^2)u_x(1.20) + field_f(ii,4) = stf_coeff * dz_f * cmplx(-om*om,0) ! (-om^2)u_z + + ! acoustic boundary point + ! note: instead of velocity as in the elastic case, in acoustic domains we would need potentials. + ! the velocity potential would be defined as: v = 1/rho grad(potential_dot) + ! thus, we would require to either change the FK formulations or to integrate velocity. + ! this will be left as a todo for future... + ! + ! as a reminder, displacement in frequency domains could be obtained by: + !if (ispec_is_acoustic(ispec)) then + ! ! displacement: u_x = -i y_1 + ! ! u_z = y_3 from Tong. (A13) + ! field_f(ii,1) = stf_coeff * dx_f * cmplx(0,-1) ! u_x = - i y_1 + ! field_f(ii,2) = stf_coeff * dz_f ! u_z = y_3 + !endif + + else + txz_f = N_mat(3,2)*coeff(1,ii) + N_mat(3,4)*coeff(2,ii) + N_mat(3,1)*C_1 ! tilde{y}_4 + tzz_f = N_mat(4,2)*coeff(1,ii) + N_mat(4,4)*coeff(2,ii) + N_mat(4,1)*C_1 ! tilde{y}_6 + + field_f(ii,1) = stf_coeff * om * ray_p * (xi1(ipt)*tzz_f - 4.0*xim(ipt)*dx_f) ! T_xx + field_f(ii,2) = stf_coeff * om * ray_p * txz_f * cmplx(0,-1) ! T_xz + field_f(ii,3) = stf_coeff * om * ray_p * tzz_f ! T_zz + endif + enddo + + ! pad negative f, and convert to time series + do ii = 2, nf2-1 + field_f(nf+2-ii,:) = conjg(field_f(ii,:)) + enddo + + field(:,:) = 0.0 + do j = 1, nvar + ! inverse FFT + call FFTinv(npow,field_f(:,j),zign_neg,dt,field(:,j),mpow) + + ! wrap around to start from t0: here one has to be careful if t0/dt is not + ! exactly an integer, assume nn > 0 + ! note: for nn == 0, nothing to wrap as field(1:npts2,j) = field(1:npts2,j) + if (nn > 0) then + dtmp(1:nn) = field(npts2-nn+1:npts2,j) + field(nn+1:npts2,j) = field(1:npts2-nn,j) + field(1:nn,j) = dtmp(1:nn) + else if (nn < 0) then + dtmp(1:nn) = field(1:nn,j) + field(1:npts-nn,j) = field(nn+1:npts,j) + field(npts-nn+1:npts,j) = dtmp(1:nn) + endif + enddo + + if (.not. compute_traction) then + !! store undersampled version of velocity FK solution + tmp_t1(:) = field(:,1)*cos(phi) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + displ(1,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(:) = field(:,1)*sin(phi) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + displ(2,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(:) = field(:,2) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + displ(3,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(:) = field(:,3)*cos(phi) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + accel(1,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(:) = field(:,3)*sin(phi) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + accel(2,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(:) = field(:,4) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + accel(3,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + else + !! compute traction + do lpts = 1, NF_FOR_STORING + sigma_rr = field(lpts,1) + sigma_rt = 0.0 + sigma_rz = field(lpts,2) + sigma_zz = field(lpts,3) + sigma_tt = bdlambdamu(ipt)*(sigma_rr+sigma_zz) + sigma_tz = 0.0 + + Txx_tmp = sigma_rr * cos(phi) * cos(phi) + sigma_tt * sin(phi) * sin(phi) + Txy_tmp = cos(phi) * sin(phi) * (sigma_rr - sigma_tt) + Txz_tmp = sigma_rz * cos(phi) + Tyy_tmp = sigma_rr * sin(phi) * sin(phi) + sigma_tt * cos(phi) * cos(phi) + Tyz_tmp = sigma_rz * sin(phi) + Tzz_tmp = sigma_zz + + !! store directly the traction + traction(1,ipt,lpts) = Txx_tmp*nmx(ipt) + Txy_tmp*nmy(ipt) + Txz_tmp*nmz(ipt) + traction(2,ipt,lpts) = Txy_tmp*nmx(ipt) + Tyy_tmp*nmy(ipt) + Tyz_tmp*nmz(ipt) + traction(3,ipt,lpts) = Txz_tmp*nmx(ipt) + Tyz_tmp*nmy(ipt) + Tzz_tmp*nmz(ipt) + enddo + + !! store undersamped version of tractions FK solution + tmp_t1(1:NF_FOR_STORING) = traction(1,ipt,1:NF_FOR_STORING) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + traction(1,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(1:NF_FOR_STORING) = traction(2,ipt,1:NF_FOR_STORING) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + traction(2,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + + tmp_t1(1:NF_FOR_STORING) = traction(3,ipt,1:NF_FOR_STORING) + call compute_spline_coef_to_store(tmp_t1, npts2, tmp_t2, tmp_c) + traction(3,ipt,1:NF_FOR_STORING) = tmp_t2(1:NF_FOR_STORING) + endif + + ! user output + if (myrank == 0 .and. npt > 1000) then + if (mod(ipt,(npt/10)) == 0) then + write(IMAIN,*) ' done ',ipt/(npt/10)*10,'% points out of ',npt + call flush_IMAIN() + endif + endif + !enddo + enddo + endif + + ! taper + if (USE_TAPERED_BEGINNING) then + ! check if taper length is reasonable compared to wavefield storage traces + if (NF_FOR_STORING > 2 * taper_nlength) then + do i = 1,taper_nlength + ! cosine taper, otherwise using a constant (1.0) instead + taper = (1.0 - cos(PI*(i-1)/taper_nlength)) * 0.5 ! between [0,1[ + if (.not. compute_traction) then + ! tapers traces + displ(:,:,i) = taper * displ(:,:,i) + accel(:,:,i) = taper * accel(:,:,i) + else + traction(:,:,i) = taper * traction(:,:,i) + endif + enddo + endif + endif + + ! free temporary arrays + deallocate(fvec,coeff, field_f, field, dtmp) + deallocate(tmp_f1, tmp_f2, tmp_f3, tmp_t1, tmp_t2, tmp_t3) + deallocate(tmp_c) + + ! user output + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) " FK computing passed " + write(IMAIN,*) + call flush_IMAIN() + endif + + end subroutine FK + +! +!------------------------------------------------------------------------------------------------- +! + + subroutine compute_N_Rayleigh(alpha, beta, mu, H, nlayer, om, ray_p, height, N_mat) + + ! assumes that height = 0 is the bottom interface + + use fk_injection, only: CUSTOM_REAL,myrank + + implicit none + + ! precision for complex + integer, parameter :: CUSTOM_CMPLX = 8 + + ! input + integer, intent(in) :: nlayer + real(kind=CUSTOM_REAL), dimension(nlayer), intent(in) :: alpha, beta, mu, H + real(kind=CUSTOM_REAL), intent(in) :: om,ray_p,height + + ! output + complex(kind=CUSTOM_CMPLX), dimension(4,4), intent(inout) :: N_mat(4,4) + + ! local vars + integer :: i, j, ilayer + complex(kind=CUSTOM_CMPLX), dimension(nlayer) :: eta_alpha, eta_beta, nu_al, nu_be + complex(kind=CUSTOM_CMPLX), dimension(4,4) :: E_mat, G_mat, P_layer + complex(kind=CUSTOM_CMPLX) :: xa, ya, xb, yb + complex(kind=CUSTOM_CMPLX) :: ca, sa, cb, sb, c1, c2 + real(kind=CUSTOM_CMPLX), dimension(nlayer) :: gamma0, gamma1 + real(kind=CUSTOM_CMPLX), dimension(nlayer) :: hh + real(kind=CUSTOM_CMPLX) :: g1, mul, two_mul, g1_sq + real(kind=CUSTOM_CMPLX) :: alphal, betal + + if (nlayer < 1) stop 'nlayer has to be larger than or equal to 1' + +! see: Tong et al. (2014), +! High-resolution seismic array imaging based on an SEM-FK hybrid method, +! GJI, 197, 369-395. +! details in the appendix + + ! note: vertical incident (p=0) is not handled in Tong et al. explanations. + ! here, it limits p to a very small value to handle the calculations + if (abs(ray_p) < 1.e-15) stop 'Invalid ray parameter p (cannot be zero) in compute_N_Rayleigh() routine' + + ! vertical wavenumbers + ! pre-computes factors + do i = 1,nlayer + ! vp and vs + alphal = alpha(i) ! vp + betal = beta(i) ! vs + + ! safety check: zero shear velocity case not handled yet + if (abs(betal) < 1.e-15) stop 'Invalid shear velocity in compute_N_Rayleigh() routine' + + ! P + ! see (A5): E_23 = -i nu_p / k = -i omega sqrt(1/alpha^2 - p^2) / k + ! factor eta_alpha = -i sqrt(1/alpha^2 - p^2) + ! = -i sqrt( (1/alpha + p) * (1/alpha - p) ) + ! = -i sqrt(1/alpha + p) * sqrt(1/alpha - p) + !org: eta_alpha(i) = -cmplx(0,1) * sqrt(1.0/alpha(i) + ray_p) * sqrt(1.0/alpha(i) - ray_p) + eta_alpha(i) = -cmplx(0,1) * sqrt( 1.0 / alphal**2 - ray_p**2 ) ! i*vertical slowness, purely imaginary + ! note: here factor nu_al = -i nu_p = -i omega sqrt(1/alpha^2 - p^2) + ! = omega (-i sqrt(1/alpha^2 - p^2) + ! = omega (eta_alpha) + nu_al(i) = om * eta_alpha(i) ! i * vertical wavenumber + + ! SV + ! see (A5): E_11 = -i nu_s / k = -i omega sqrt(1/beta^2 - p^2) / k + ! factor eta_beta = -i sqrt(1/beta^2 - p^2) + ! = -i sqrt( (1/beta + p) * (1/beta - p) ) + ! = -i sqrt(1/beta + p) * sqrt(1/beta - p) + !org: eta_beta(i) = -cmplx(0,1) * sqrt(1.0/beta(i) + ray_p) * sqrt(1.0/beta(i) - ray_p) + eta_beta(i) = -cmplx(0,1) * sqrt( 1.0 / betal**2 - ray_p**2 ) + + ! note: here factor nu_be = -i nu_s = -i omega sqrt(1/beta^2 - p^2) + ! = omega (-i sqrt(1/beta^2 - p^2) + ! = omega (eta_beta) + nu_be(i) = om * eta_beta(i) ! i * vertical wavenumber + + ! auxiliary variables + gamma0(i) = 2.0 * betal**2 * ray_p**2 + gamma1(i) = 1.0 - 1.0/gamma0(i) + enddo + + ! initializes matrix + E_mat(:,:) = (0.0,0.0) + + ! Tong et al. (2014), appendix (A10) E_0: + ! note: E_mat is not omega dependent + ! + ! (A5): E_11 = -i nu_s / k + ! = -i omega * sqrt(1/beta^2 - p^2) / ( p * omega) ,with p = k/omega -> k = p * omega + ! = -i sqrt(1/beta^2 - p^2) / p + ! = eta_beta / p + E_mat(1,1) = eta_beta(nlayer) / ray_p + E_mat(1,2) = -E_mat(1,1) + E_mat(1,3) = (1.0,0.0) + E_mat(1,4) = (1.0,0.0) + + ! (A5): E_23 = -i nu_p / k + ! = -i omega * sqrt(1/alpha^2 - p^2) / (p * omega) ,with p = k/omega -> k = p * omega + ! = -i sqrt(1/alpha^2 - p^2) / p + ! = eta_alpha / p + E_mat(2,1) = (1.0,0.0) + E_mat(2,2) = (1.0,0.0) + E_mat(2,3) = eta_alpha(nlayer) / ray_p + E_mat(2,4) = -E_mat(2,3) + + ! pre-computed factor for half-space (layer with index nlayer) + mul = mu(nlayer) + ! 2 mu + two_mul = 2.0 * mul + + ! note: wavenumber k is defined as p = k / omega -> k = p * omega + ! for frequency omega==0, wavenumber k becomes zero and thus a factor 1/k becomes undefined. + ! to avoid such behavior, the factor k in the expressions for E below will be cancelled out. + ! together with the multiplication by propagation matrix P, where elements can have a factor k, this should be ok. + ! unfortunately, the original paper by Tong et al. doesn't mention such a case and assumes k is non-zero. + + ! (A5): E_31 = 2 k mu gamma_1 + ! = 2 p omega mu gamma_1 ,with wavenumber k: from ray p = k / omega -> k = p * omega + !org: E_mat(3,1) = 2.0 * mu(nlayer) * gamma1(nlayer) ! takes out factor k + E_mat(3,1) = two_mul * gamma1(nlayer) + E_mat(3,2) = E_mat(3,1) + + ! (A5): E_33 = - 2 i mu nu_p + ! = 2 mu (-i) omega sqrt(1/alpha^2 - p^2) + ! = 2 mu omega (-i sqrt(1/alpha^2 - p^2)) + ! = 2 mu omega (eta_alpha) + ! or + ! = 2 mu (eta_alpha) k / p ,with p = k/omega -> omega = k / p + !org: E_mat(3,3) = 2.0 * mu(nlayer) * eta_alpha(nlayer) / ray_p ! takes out factor k + E_mat(3,3) = two_mul * eta_alpha(nlayer) / ray_p + E_mat(3,4) = -E_mat(3,3) + + ! (A5): E_41 = - 2 i mu nu_s + ! = 2 mu (-i) omega sqrt(1/beta^2 - p^2) + ! = 2 mu omega (-i sqrt(1/beta^2 - p^2)) + ! = 2 mu omega (eta_beta) + ! or + ! = 2 mu (eta_beta) k / p ,with p = k/omega -> omega = k / p + !org: E_mat(4,1) = 2.0 * mu(nlayer) * eta_beta(nlayer) / ray_p ! takes out factor k + E_mat(4,1) = two_mul * eta_beta(nlayer) / ray_p + E_mat(4,2) = -E_mat(4,1) + + ! (A5): E_43 = 2 k mu gamma_1 + ! = E_31 = E_32 + E_mat(4,3) = E_mat(3,1) + E_mat(4,4) = E_mat(3,1) + + if (height > sum(H(1:nlayer-1))) then + print *,'FK error: rank ',myrank + print *,' Z point is located in the air above the surface rather than in the solid!' + print *,' current z :', height, ' max z allowed : ', sum(H(1:nlayer-1)) + stop 'FK invalid height' + endif + + ! figure out the location z with respect to layer stack + if (height <= 0.0) then + ! in lower half space + G_mat(:,:) = (0.0,0.0) + ! incident, up-going S-wave + ! (A5): Gamma_11 = e^(-i nu_s z) = e^( (-i nu_s) * z ) = e^(nu_be * z) + G_mat(1,1) = exp(nu_be(nlayer) * height) + ! reflected, down-going S-wave + ! (A5): Gamma_22 = e^(i nu_s z) = e^( -(-i nu_s) * z ) = e^(-nu_be * z) + G_mat(2,2) = exp(-nu_be(nlayer) * height) + ! incident, up-going P-wave + ! (A5): Gamma_33 = e^(-i nu_p z) = e^( (-i nu_p) * z ) = e^(nu_al * z) + G_mat(3,3) = exp(nu_al(nlayer) * height) + ! reflected, down-going P-wave + ! (A5): Gamma_44 = e^(i nu_p z) = e^( -(-i nu_p) * z ) = e^(-nu_al * z) + G_mat(4,4) = exp(-nu_al(nlayer) * height) + + ! resulting matrix + N_mat = matmul(E_mat,G_mat) + + else + ! in layers + ! determines layer in which the point (given by height) lies + ! note: indexing assumes that last layer (nlayer) being the bottom, lower halfspace, + ! and the first layer (1) being at the top surface + hh(:) = H(:) + ilayer = nlayer + do j = nlayer-1 , 1 , -1 + if (height <= sum(H(j:nlayer-1))) then + ilayer = j; exit + endif + enddo + + ! updates point's layer thicknesses + hh(ilayer+1:nlayer-1) = H(ilayer+1:nlayer-1) + hh(ilayer) = height - sum(H(ilayer+1:nlayer-1)) + + if (hh(ilayer) < 0.0) then + print *,'Error: rank ',myrank,' has invalid point height ',hh(ilayer),' at layer ',ilayer + stop 'Error setting layer thickness' + endif + + !debug + !print *,'debug: height ',height,'layer ',ilayer,nlayer,'H',H(:),'sum',sum(H(ilayer:nlayer-1)),'h',hh(:) + + ! compute propagation matrices + N_mat(:,:) = E_mat(:,:) + + do j = nlayer-1, ilayer, -1 + ! matrix variables + ! C_alpha = cos(nu_p h) = cos( omega sqrt(1/alpha^2 - p^2) h ) + ! S_alpha = -sin(nu_p h) = -sin( omega sqrt(1/alpha^2 - p^2) h ) + ! with nu_p h = omega sqrt(1/alpha^2 - p^2) h + ! + ! note: there might be some sign conflict and confusion between sin and sinh in the definition after (A9) of Tong et al. + ! instead of S_alpha = -sin(nu_p h) as stated in the paper, here sa becomes [i sin(nu_p h)] as imaginary number. + ! + c1 = nu_al(j) * hh(j) + ! cos(nu_p h) = [ e^(i nu_p h) + e^(-i nu_p h) ] / 2 + ! = [ e^(-nu_al * h) + e^(nu_al * h) ] / 2 + ca = (exp(c1) + exp(-c1))/2.0 + ! sin(nu_p h) = [ e^(i nu_p h) - e^(-i nu_p h) ] / 2i + ! = [ e^(-nu_al * h) - e^(nu_al * h) ] / 2i + ! = i [ e^(nu_al * h) - e^(-nu_al * h) ] / 2 + sa = (exp(c1) - exp(-c1))/2.0 ! imaginary part + + ! X_alpha = - i nu_p S_alpha / ( omega p) + ! = -i omega sqrt(1/alpha^2 - p^2) S_alpha / (omega p ) + ! = [ -i sqrt(1/alpha^2 - p^2) ] S_alpha / p + ! = eta_alpha S_alpha / p + xa = eta_alpha(j) * sa / ray_p + + ! Y_alpha = i omega p S_alpha / nu_p + ! = p S_alpha i omega / ( omega sqrt(1/alpha^2 - p^2) ) + ! = p S_alpha i / ( sqrt(1/alpha^2 - p^2) ) + ! = p S_alpha 1 / (-i sqrt(1/alpha^2 - p^2) ) + ! = p S_alpha / eta_alpha + ya = ray_p * sa / eta_alpha(j) + + ! C_beta = cos(nu_s h) = cos( omega sqrt(1/beta^2 - p^2) h ) + ! S_beta = -sin(nu_s h) = -sin( omega sqrt(1/beta^2 - p^2) h ) + ! with nu_s h = omega sqrt(1/beta^2 - p^2) h + ! + ! note: for sb, see same remark as above for sa + c2 = nu_be(j) * hh(j) + cb = (exp(c2) + exp(-c2))/2.0 ! cos(nu_s h) + sb = (exp(c2) - exp(-c2))/2.0 ! sin(nu_s h) imaginary part + + ! X_beta = - i nu_s S_beta / ( omega p) + ! Y_beta = i omega p S_beta / nu_s + xb = eta_beta(j) * sb / ray_p + yb = ray_p * sb / eta_beta(j) + + ! layer factors + g1 = gamma1(j) + mul = mu(j) + ! pre-computed factors + ! gamma1^2 + g1_sq = g1 * g1 + ! 2 mu + two_mul = 2.0 * mul ! note: leaving out factor k, only 2 mu + + ! Tong et al. (2014), appendix (A8) + ! propagation matrix P_n + P_layer(1,1) = ca - g1*cb + P_layer(1,2) = xb - g1*ya + !org: P_layer(1,3) = (ya - xb)/(2*mul) ! misses factor 1/k + ! P_layer(1,4) = (cb - ca)/(2*mul) ! misses factor 1/k + P_layer(1,3) = (ya - xb) / two_mul + P_layer(1,4) = (cb - ca) / two_mul + + P_layer(2,1) = xa - g1*yb + P_layer(2,2) = cb - g1*ca + !org: P_layer(2,3) = (ca - cb)/(2*mul) ! misses factor 1/k + ! P_layer(2,4) = (yb - xa)/(2*mul) ! misses factor 1/k + P_layer(2,3) = (ca - cb) / two_mul + P_layer(2,4) = (yb - xa) / two_mul + + !org: P_layer(3,1) = 2*mul * (xa - g1**2 * yb) ! misses factor k + ! P_layer(3,2) = 2*mul * g1 * (cb - ca) ! misses factor k + P_layer(3,1) = two_mul * (xa - g1_sq * yb) + P_layer(3,2) = two_mul * g1 * (cb - ca) + + P_layer(3,3) = ca - g1*cb + P_layer(3,4) = g1*yb - xa + + !org: P_layer(4,1) = 2*mul * g1 * (ca - cb) ! misses factor k + ! P_layer(4,2) = 2*mul * (xb - g1**2 * ya) ! misses factor k + P_layer(4,1) = two_mul * g1 * (ca - cb) + P_layer(4,2) = two_mul * (xb - g1_sq * ya) + P_layer(4,3) = g1*ya - xb + P_layer(4,4) = cb - g1*ca + + !debug + !if (myrank == 0) print *,'debug: j,g1,g1_sq,mul,two_mul,xa,xb,ya,yb,ca,cb,ray_p,om', & + ! j,g1,g1_sq,mul,two_mul,xa,xb,ya,yb,ca,cb,ray_p,om + !if (myrank == 0) print *,'debug: j,P_layer',j,P_layer(:,:) + + ! resulting matrix + N_mat = gamma0(j) * matmul(P_layer,N_mat) + enddo + endif + + ! debug + !if (myrank == 0) then + ! print *,'debug: N_mat ' + ! do j = 1,4 + ! print *,N_mat(:,j) + ! enddo + !endif + !stop + + end subroutine compute_N_Rayleigh + +! +!------------------------------------------------------------------------------------------------- +! + + subroutine FFT(npow,xi,zign,dtt,mpow) + +! Fourier transform +! This inputs AND outputs a complex function. +! The convention is FFT --> e^(-iwt) +! numerical factor for Plancherel theorem: planch_fac = dble(NPT * dt * dt) + + use fk_injection, only: CUSTOM_REAL + + implicit none + + integer, parameter :: CUSTOM_CMPLX = 8 + + integer,intent(in) :: npow + + real(kind=CUSTOM_REAL),intent(in) :: zign,dtt + + complex(kind=CUSTOM_CMPLX),dimension(*),intent(inout) :: xi + +!! DK DK here is the hardwired maximum size of the array +!! DK DK Aug 2016: if this routine is called many times (for different mesh points at which the SEM is coupled with FK) +!! DK DK Aug 2016: this should be moved to the calling program and precomputed once and for all + real(kind=CUSTOM_REAL),intent(in) :: mpow(30) + + ! local parameters + integer :: lblock,k,FK,jh,ii,istart + integer :: l,iblock,nblock,i,lbhalf,j,lx + complex(kind=CUSTOM_CMPLX) :: wk, hold, q + real(kind=CUSTOM_REAL) :: flx,inv_of_flx,v + ! PI = acos(-1.d0) + real(kind=CUSTOM_REAL), parameter :: TWO_PI = 2.0_CUSTOM_REAL * acos(-1.d0) + +!! DK DK added this sanity check + if (npow > 30) stop 'error: the FK FTT routine has an hardwired maximum of 30 levels' +!! DK DK in any case the line below imposes a maximum of 31, otherwise the integer 2**n will overflow + + lx = 2**npow + + do l = 1,npow + + nblock = 2**(l-1) + lblock = lx/nblock + lbhalf = lblock/2 + + k = 0 + + do iblock = 1,nblock + + FK = k + flx = lx + + v = zign * TWO_PI * FK / flx ! Fourier convention + + ! - sign: MATLAB convention: forward e^{-i om t} + ! + sign: engineering convention: forward e^{i om t} + wk = cmplx(cos(v),-sin(v)) ! sign change to -sin(v) or sin(v) + istart = lblock*(iblock-1) + + do i = 1,lbhalf + j = istart+i + jh = j+lbhalf + q = xi(jh)*wk + xi(jh) = xi(j)-q + xi(j) = xi(j)+q + enddo + + do i = 2,npow + ii = i + if (k < mpow(i)) goto 4 + k = k-mpow(i) + enddo + + 4 k = k+mpow(ii) + + enddo + enddo + + k = 0 + + do j = 1,lx + if (k < j) goto 5 + + hold = xi(j) + xi(j) = xi(k+1) + xi(k+1) = hold + +5 do i = 1,npow + ii = i + if (k < mpow(i)) goto 7 + k = k-mpow(i) + enddo + +7 k = k+mpow(ii) + enddo + + ! final steps deal with dt factors + if (zign > 0.) then + ! FORWARD FFT + xi(1:lx) = xi(1:lx) * dtt ! multiplication by dt + + else + ! REVERSE FFT + flx = flx*dtt + inv_of_flx = 1._CUSTOM_REAL / flx + +!! DK DK Aug 2016: changed to multiplication by the precomputed inverse to make the routine faster +! xi(1:lx) = xi(1:lx) / flx ! division by dt + xi(1:lx) = xi(1:lx) * inv_of_flx ! division by dt + + endif + + end subroutine FFT + +! +!------------------------------------------------------------------------------------------------- +! + + subroutine FFTinv(npow,s,zign,dtt,r,mpow) + +! inverse Fourier transform -- calls FFT + + use fk_injection, only: CUSTOM_REAL + + implicit none + + integer, parameter :: CUSTOM_CMPLX = 8 + + integer,intent(in) :: npow + real(kind=CUSTOM_REAL),intent(in) :: dtt,zign + + complex(kind=CUSTOM_CMPLX), intent(inout) :: s(*) + real(kind=CUSTOM_REAL), intent(out) :: r(*) ! note that this is real, not double precision + +!! DK DK here is the hardwired maximum size of the array +!! DK DK Aug 2016: if this routine is called many times (for different mesh points at which the SEM is coupled with FK) +!! DK DK Aug 2016: this should be moved to the calling program and precomputed once and for all + real(kind=CUSTOM_REAL),intent(in) :: mpow(30) + + ! local parameters + integer :: nsmp, nhalf + + nsmp = 2**npow + nhalf = nsmp/2 + + call rspec(s,nhalf) ! restructuring + call FFT(npow,s,zign,dtt,mpow) ! Fourier transform + + r(1:nsmp) = real(s(1:nsmp)) ! take the real part + + end subroutine FFTinv + +! +!------------------------------------------------------------------------------------------------- +! + + subroutine rspec(s,np2) + + implicit none + + integer, parameter :: CUSTOM_CMPLX = 8 + + complex(kind=CUSTOM_CMPLX),intent(inout) :: s(*) + integer, intent(in) :: np2 + + ! local parameters + integer :: n,n1,i + + n = 2*np2 + n1 = np2+1 + + s(n1) = 0.0 + s(1) = cmplx(real(s(1)),0.0) + + do i = 1,np2 + s(np2+i) = conjg(s(np2+2-i)) + enddo + + end subroutine rspec + +! +!------------------------------------------------------------------------------------------------- +! + + subroutine find_size_of_working_arrays(deltat, freq_sampling_fk, tmax_fk, NF_FOR_STORING, & + NF_FOR_FFT, NPOW_FOR_INTERP, NP_RESAMP, DF_FK) + + use fk_injection, only: CUSTOM_REAL, NSTEP + + implicit none + real(kind=CUSTOM_REAL),intent(in) :: deltat,freq_sampling_fk + real(kind=CUSTOM_REAL),intent(inout) :: tmax_fk + real(kind=CUSTOM_REAL),intent(inout) :: DF_FK + integer, intent(inout) :: NF_FOR_STORING, NF_FOR_FFT, NPOW_FOR_INTERP, NP_RESAMP + ! local parameter + real(kind=CUSTOM_REAL) :: df, dt_min_fk, Frq_ech_Fk + real(kind=CUSTOM_REAL) :: tmax_samp, tmax_use + + ! sampling frequency to store fk solution + Frq_ech_Fk = freq_sampling_fk + + ! sampling time step to store fk solution + dt_min_fk = 1.0 / Frq_ech_Fk + + ! compute resampling rate + NP_RESAMP = floor(dt_min_fk / deltat) + + ! checks sampling rate + if (NP_RESAMP == 0) NP_RESAMP = 1 + + ! update dt for fk with respect to integer np_resampling + dt_min_fk = NP_RESAMP * deltat !! this is the time step sampling for FK storage + + !! time window lengths + ! maximum simulation time length would be: + ! tmax_simulation = NSTEP * deltat + ! FK signal trace is symmetric after inverse FFT, thus NF_FOR_STORING must match at least match NSTEP * 2 with resampling rate. + ! maximum time length with respect to simulation NSTEP and resampling rate: + tmax_samp = (2*NSTEP/NP_RESAMP + 1) * dt_min_fk + + ! take user defined window length as initial default + tmax_use = tmax_fk + + ! limits window length to actual simulation length + if (tmax_use > tmax_samp) tmax_use = tmax_samp + + ! compute number of time steps to store + NF_FOR_STORING = ceiling( tmax_use / dt_min_fk) + + ! in power of two + NF_FOR_STORING = ceiling(log(real(NF_FOR_STORING))/log(2.)) + + ! multiply by 2 in order to do an inverse FFT + NF_FOR_FFT = 2**(NF_FOR_STORING+1) + + NPOW_FOR_INTERP = NF_FOR_STORING+1 + NF_FOR_STORING = 2**NF_FOR_STORING + + ! now we have this new time window + tmax_fk = dt_min_fk * (NF_FOR_FFT - 1) + + ! step in frequency for fk + df = 1.0 / tmax_fk + DF_FK = df + + !debug + !print *,'debug: tmax ',tmax_fk,tmax_samp,tmax_use,'np_resamp',NP_RESAMP,'NSTEP',NSTEP,2*NSTEP/NP_RESAMP, & + ! 'NF',NF_FOR_STORING,ceiling( tmax_use / dt_min_fk) + + end subroutine find_size_of_working_arrays + + +! +!------------------------------------------------------------------------------------------------- +! + +!! ################# INTERPOLATION ROUTINES IN TIME DOMAIN ###################################### + +!! compute and store spline coefficients + + subroutine compute_spline_coef_to_store(Sig, npts, spline_coeff, tmp_c) + + use fk_injection, only: CUSTOM_REAL + + implicit none + + integer, intent(in) :: npts + real(kind=CUSTOM_REAL), dimension(npts), intent(in) :: Sig + real(kind=CUSTOM_REAL), dimension(npts), intent(inout) :: spline_coeff + double precision, dimension(npts),intent(out) :: tmp_c ! temporary work array + + !! computation in double precision + double precision,parameter :: error = 1.d-24 + double precision,parameter :: z1 = dsqrt(3.d0) - 2.d0 + double precision :: zn, sumc + integer :: i, n_init + + ! Compute pole value + tmp_c(:) = dble(Sig(:)) * (1.d0 - z1) * ( 1.d0 - 1.d0/z1) + + ! Initialisation causal filter + n_init = ceiling(log(error)/log(abs(z1))) + + ! check limits: by default is n_init==42, for very short test simulations this might be larger than npts + if (n_init > npts) n_init = npts + + sumc = tmp_c(1) + zn = z1 + do i = 1,n_init + sumc = sumc + zn * tmp_c(i) + zn = zn * z1 + enddo + tmp_c(1) = sumc + + ! Causal filter + do i = 2,npts + tmp_c(i) = tmp_c(i) + z1 * tmp_c(i-1) + enddo + + ! Initialisation anti-causal filter + tmp_c(npts) = ( z1 / (z1 - 1.d0) ) * tmp_c(npts) + do i = npts-1,1,-1 + tmp_c(i) = z1 * (tmp_c(i+1) - tmp_c(i)) + enddo + + !! store spline coeff in CUSTOM_REAL precision + spline_coeff(:) = tmp_c(:) + + end subroutine compute_spline_coef_to_store + + +! +!------------------------------------------------------------------------------------------------- +! + +!! VM VM READING INPUT FILE FOR FK MODEL + + subroutine ReadFKModelInput() + + !use specfem_par + !use specfem_par_coupling + use fk_injection, only: CUSTOM_REAL, TINYVAL, alpha_FK, beta_FK, & + rho_FK, mu_FK, h_FK, phi_FK, theta_FK, & + nlayer, IMAIN, myrank, & + xx0, yy0, zz0, tt0, ff0, freq_sampling_fk, & + tmax_fk, amplitude_fk, Z_REF_for_FK, NSTEP, deltat, & + type_kpsv_fk, FKMODEL_FILE + + implicit none + + !real(kind=CUSTOM_REAL), intent(in) :: Xmin_box, Xmax_box, Ymin_box, Ymax_box, Zmin_box, Zmax_box + integer :: ioerr + character(len=100) :: keyword, keyvalue, line + character(len=100) :: keyword_tmp, incident_wave + + real(kind=CUSTOM_REAL) :: rho_layer, vp_layer, vs_layer, ztop_layer + real(kind=CUSTOM_REAL) :: Radius_box, wave_length_at_bottom + real(kind=CUSTOM_REAL), dimension(:), allocatable :: rho_fk_input, vp_fk_input, vs_fk_input, ztop_fk_input + integer, dimension(:), allocatable :: ilayer_fk_input + integer :: ilayer,ier + !logical :: position_of_wavefront_not_read + + !!-------------------------------------------------------------- + ! # model description : + ! NLAYER n # number of layers + ! LAYER 1 rho, vp ,vs, ztop + ! LAYER 2 rho, vp, vs, ztop + ! ... + ! LAYER n rho, vp, vs, ztop # homogenoeus half-space + ! # + ! # incident wave description: + ! INCIDENT_WAVE "p" or "sv" + ! BACK_AZITUTH bazi + ! INCIDENCE inc + ! ORIGIN_WAVEFRONT xx0, yy0, zz0 + ! ORIGIN_TIME tt0 + ! FREQUENCY_MAX ff0 + ! FREQUENCY_SAMPLING freq_sampling_fk + ! TIME_WINDOW tmax_fk + ! AMPLITUDE amplitude_fk + !!---------------------------------------------------------------- + + !! set default values + tt0 = 0.0 ! origin time + tmax_fk = 128.0 ! time window length + + ff0 = 0.1 ! frequency maximum in Hz + freq_sampling_fk = 10.0 ! frequency sampling rate in Hz + amplitude_fk = 1.0 ! amplitude factor + + type_kpsv_fk = 1 ! 1 == P-wave / 2 == SV-wave + + !position_of_wavefront_not_read = .true. + + !! READING input file + open(85,file=trim(FKMODEL_FILE)) + do + read(85, fmt='(a100)',iostat=ioerr) line + !call remove_begin_blanks(line) + if (ioerr < 0) exit + if (len(trim(line)) < 1 .or. line(1:1) == '#') cycle + read(line,*) keyword, keyvalue + + select case(trim(keyword)) + case('NLAYER') + read(line, *) keyword_tmp, nlayer + + allocate(alpha_FK(nlayer), beta_FK(nlayer), rho_FK(nlayer), mu_FK(nlayer), h_FK(nlayer),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating arrays 2226') + alpha_FK(:) = 0._CUSTOM_REAL; beta_FK(:) = 0._CUSTOM_REAL; rho_FK(:) = 0._CUSTOM_REAL + mu_FK(:) = 0._CUSTOM_REAL; h_FK(:) = 0._CUSTOM_REAL + + allocate(rho_fk_input(nlayer), & + vp_fk_input(nlayer), & + vs_fk_input(nlayer), & + ztop_fk_input(nlayer+1), & + ilayer_fk_input(nlayer+1),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating arrays 2227') + rho_fk_input(:) = 0._CUSTOM_REAL + vp_fk_input(:) = 0._CUSTOM_REAL + vs_fk_input(:) = 0._CUSTOM_REAL + ztop_fk_input(:) = 0._CUSTOM_REAL + ilayer_fk_input(:) = -1 + + case('LAYER') + read(line, *) keyword_tmp, ilayer, rho_layer, vp_layer, vs_layer, ztop_layer + + ilayer_fk_input(ilayer) = ilayer + rho_fk_input(ilayer) = rho_layer + vp_fk_input(ilayer) = vp_layer + vs_fk_input(ilayer) = vs_layer + ztop_fk_input(ilayer) = ztop_layer + + ! acoustic case: FK for zero-shear velocity not handled, here we set it to a (very) small value + if (abs(vs_fk_input(ilayer)) < TINYVAL) vs_fk_input(ilayer) = 1.e-5 + + case('INCIDENT_WAVE') + read(line,*) keyword_tmp, incident_wave + + select case(trim(incident_wave)) + case ('p', 'P') + type_kpsv_fk = 1 + case('sv','SV') + type_kpsv_fk = 2 + case default + type_kpsv_fk = 1 + end select + + case('BACK_AZIMUTH') + read(line,*) keyword_tmp, phi_FK + phi_FK = - phi_FK - 90.0 + + case('AZIMUTH') + read(line,*) keyword_tmp, phi_FK + phi_FK = 90.0 - phi_FK + + case('TAKE_OFF') + read(line,*) keyword_tmp, theta_FK + + case('ORIGIN_WAVEFRONT') + read(line,*) keyword_tmp, xx0, yy0, zz0 + !position_of_wavefront_not_read = .false. + + case('ORIGIN_TIME') + read(line,*) keyword_tmp, tt0 + + case('NSTEP') + read(line,*) keyword_tmp, NSTEP + + case('deltat') + read(line,*) keyword_tmp, deltat + + case('FREQUENCY_MAX') + read(line,*) keyword_tmp, ff0 + + case('FREQUENCY_SAMPLING') + read(line,*) keyword_tmp, freq_sampling_fk + + case('TIME_WINDOW') + read(line,*) keyword_tmp, tmax_fk + + case('AMPLITUDE') + read(line,*) keyword_tmp, amplitude_fk + + end select + !!------------------------------------------------------------------------------------------------------ + enddo + + if (allocated(ilayer_fk_input)) then + + ilayer_fk_input(nlayer+1) = ilayer_fk_input(nlayer) + ztop_fk_input(nlayer+1) = ztop_fk_input(nlayer) + + ! Z coordinate reference for top of lower-halfspace (last layer input) + Z_REF_for_FK = ztop_fk_input(nlayer) + + do ilayer = 1, nlayer + alpha_FK(ilayer) = vp_fk_input(ilayer) + beta_FK(ilayer) = vs_fk_input(ilayer) + rho_FK(ilayer) = rho_fk_input(ilayer) + + mu_FK(ilayer) = rho_fk_input(ilayer) * vs_fk_input(ilayer)**2 + h_FK(ilayer) = ztop_fk_input(ilayer) - ztop_fk_input(ilayer+1) + + ! checks if layer has values + if (ilayer_fk_input(ilayer) == -1) then + write(*,*) " ERROR READING FK INPUT FILE " + write(*,*) " MISSING LAYER ", ilayer + stop 'Missing layer in FK model' + endif + + ! checks vp,rho are strictly positive + if (alpha_FK(ilayer) <= 0.0_CUSTOM_REAL .or. rho_FK(ilayer) <= 0.0_CUSTOM_REAL) then + print *,'Error: invalid elastic material for FK coupling in layer ',ilayer + print *,' vp = ',alpha_FK(ilayer),' vs = ',beta_FK(ilayer),' rho = ',rho_FK(ilayer) + print *,'vp & rho must be strictly positive' + stop 'Invalid material for FK coupling found in ReadFKModelInput() routine' + endif + + ! checks if vs is strictly positive for SV coupling + if (type_kpsv_fk == 2) then + if (beta_FK(ilayer) <= 0.0_CUSTOM_REAL) then + print *,'Error: invalid elastic material for FK coupling in layer ',ilayer + print *,' vp = ',alpha_FK(ilayer),' vs = ',beta_FK(ilayer),' rho = ',rho_FK(ilayer) + print *,'vp, vs & rho must be strictly positive (in particular vs must be > 0) for SV incident wave' + stop 'Invalid elastic material for FK coupling found in ReadFKModelInput() routine' + endif + endif + + enddo + + deallocate(ilayer_fk_input, rho_fk_input, vp_fk_input, vs_fk_input, ztop_fk_input) + + else + + write(*,*) " ERROR READING FK INPUT FILE " + write(*,*) " NOT BE ABLE TO READ MODEL PROPERTIES " + stop 'Error reading FK input file' + + endif + + !! compute position of wave front + !if (position_of_wavefront_not_read) then + ! ! sets center point of box + ! xx0 = 0.5*(Xmin_box + Xmax_box) + ! yy0 = 0.5*(Ymin_box + Ymax_box) + + ! Radius_box = sqrt( (Xmin_box - xx0)**2 + (Ymin_box - yy0)**2) + + ! if (type_kpsv_fk == 1) then + ! ! P-wave + ! wave_length_at_bottom = alpha_FK(nlayer) / ff0 ! vp + ! else if (type_kpsv_fk == 2) then + ! ! SV-wave + ! wave_length_at_bottom = beta_FK(nlayer) / ff0 ! vs + ! endif + + ! ! depth position below bottom + ! zz0 = Zmin_box - Radius_box * sin ( abs(theta_FK) * (PI / 180.d0) ) & + ! - 3.0 * wave_length_at_bottom * cos ( abs(theta_FK) * (PI / 180.d0) ) + !endif + + ! user output + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) "**********************************************" + write(IMAIN,*) " USING FK INJECTION TECHNIQUE " + write(IMAIN,*) "**********************************************" + write(IMAIN,*) + write(IMAIN,*) " Model : " , nlayer , " layers " + write(IMAIN,*) + do ilayer = 1, nlayer + write(IMAIN,'(a9, i3, 3(a6,2x,f8.3), 3x, a9, f18.5 )') & + " layer ", ilayer, & + " rho =", rho_FK(ilayer), & + " vp =", alpha_FK(ilayer), & + " vs =", beta_FK(ilayer), & + " Height =", h_FK(ilayer) + enddo + write(IMAIN,*) + write(IMAIN,*) " FK (azimuth) angle phi = ", phi_FK,'(deg)' + write(IMAIN,*) " FK take-off angle theta = ", theta_FK,'(deg)' + write(IMAIN,*) + write(IMAIN,*) " Origin wavefront point FK : ", xx0, yy0, zz0 + write(IMAIN,*) " Time shift FK : ", tt0 + write(IMAIN,*) " Time discretization in SEM : ", deltat + write(IMAIN,*) " Number of steps in SEM : ", NSTEP + write(IMAIN,*) " Z reference for FK routine : ", Z_REF_for_FK + write(IMAIN,*) + write(IMAIN,*) " Time window for FK computing : ", tmax_fk + write(IMAIN,*) " Frequency max : ", ff0 + write(IMAIN,*) " Frequency sampling rate : ", freq_sampling_fk + write(IMAIN,*) " Amplitude factor : ", amplitude_fk + write(IMAIN,*) + !write(IMAIN,*) " Model boundary box : Xmin/Xmax = ",Xmin_box,Xmax_box + !write(IMAIN,*) " : Ymin/Ymax = ",Ymin_box,Ymax_box + !write(IMAIN,*) " : Zmin/Zmax = ",Zmin_box,Zmax_box + write(IMAIN,*) + write(IMAIN,*) " Type of incoming wave (1=P), (2=SV) : ",type_kpsv_fk + write(IMAIN,*) + call flush_IMAIN() + endif + + end subroutine ReadFKModelInput diff --git a/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/make_fk_injection.sh b/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/make_fk_injection.sh new file mode 100644 index 000000000..e9dd1b09d --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/make_fk_injection.sh @@ -0,0 +1,6 @@ +MPIFC=mpif90 +FLAGS='-g -xHost -O3' + +${MPIFC} ${FLAGS} -o compute_fk_injection_field coupling_fk.f90 utils.f90 compute_fk_injection_field.f90 + +${MPIFC} ${FLAGS} -o compute_fk_receiver coupling_fk.f90 utils.f90 compute_fk_receiver.f90 diff --git a/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/utils.f90 b/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/utils.f90 new file mode 100644 index 000000000..e15e56c29 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/fk_coupling/utils.f90 @@ -0,0 +1,72 @@ + subroutine flush_IMAIN() + + use fk_injection, only: IMAIN + + implicit none + + ! only main process writes out to main output file + ! file I/O in Fortran is buffered by default + ! + ! note: Fortran2003 includes a FLUSH statement + ! which is implemented by most compilers by now + ! + ! otherwise: + ! a) comment out the line below + ! b) try to use instead: call flush(IMAIN) + + flush(IMAIN) + + end subroutine flush_IMAIN + +! version without rank number printed in the error message + + subroutine exit_MPI_without_rank(error_msg) + + use mpi + !use fk_injection + + implicit none + + character(len=*) :: error_msg + integer :: ier + + ! write error message to screen + write(*,*) error_msg(1:len(error_msg)) + write(*,*) 'Error detected, aborting MPI...' + + ! flushes possible left-overs from print-statements + call flush_stdout() + + ! abort execution + call MPI_ABORT(MPI_COMM_WORLD,30,ier) + stop 'error, program ended in abort_mpi' + + end subroutine exit_MPI_without_rank + + subroutine flush_stdout() + +! flushes possible left-overs from print-statements + + implicit none + + logical :: is_connected + + ! note: Cray systems don't flush print statements before ending with an MPI + ! abort, + ! which often omits debugging statements with print before it. + ! + ! to check which unit is used for standard output, one might also use a + ! Fortran2003 module iso_Fortran_env: + ! use, intrinsic :: iso_Fortran_env, only: output_unit + + ! checks default stdout unit 6 + inquire(unit=6,opened=is_connected) + if (is_connected) & + flush(6) + + ! checks Cray stdout unit 101 + inquire(unit=101,opened=is_connected) + if (is_connected) & + flush(101) + + end subroutine flush_stdout diff --git a/EXAMPLES/applications/wavefield_discontinuity/plot/plot_utils_injection.py b/EXAMPLES/applications/wavefield_discontinuity/plot/plot_utils_injection.py new file mode 100644 index 000000000..3a0fe0dc1 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/plot/plot_utils_injection.py @@ -0,0 +1,255 @@ +import numpy as np +import sys, os +import matplotlib.pyplot as plt +from scipy import signal +from scipy.interpolate import interp1d + +def split_line(line): + # split line with continuous multiple spaces + # return a list of strings + return [_ for _ in line.strip().split(' ') if _ != ''] + +class WaveformSection: + def __init__(self): + self.waveforms = None + self.t = None + self.coord_list = None + + def __sub__(self, other): + sub = other.interpolate(self.t) + for i_wave in range(len(sub.waveforms)): + sub.waveforms[i_wave] = self.waveforms[i_wave] - sub.waveforms[i_wave] + return sub + + def shift(self, dt): + self.t += dt + + def scale(self, fac): + for i_wave in range(len(self.waveforms)): + self.waveforms[i_wave] = self.waveforms[i_wave] * fac + + def difference_ratio(self, other, time_range=None, offset=None): + if offset is None: offset = [0.0] * len(self.waveforms) + if time_range is None: time_range = (self.t[0], self.t[-1]) + ratio_list = [] + sub = other.interpolate(self.t) + for i_wave in range(len(sub.waveforms)): + ind = np.logical_and((self.t-offset[i_wave] >= time_range[0]), (self.t-offset[i_wave] <= time_range[1])) + diff = self.waveforms[i_wave][ind] - sub.waveforms[i_wave][ind] + ratio = np.sum(diff * diff) / ( + np.sqrt(np.sum(self.waveforms[i_wave][ind] * self.waveforms[i_wave][ind]))* + np.sqrt(np.sum(sub.waveforms[i_wave][ind] * sub.waveforms[i_wave][ind]))) + ratio_list.append(ratio) + return ratio_list + + def time_shift(self, other, time_range=None, offset=None): + if offset is None: offset = [0.0] * len(self.waveforms) + if time_range is None: time_range = (self.t[0], self.t[-1]) + shift_list = [] + sub = other.interpolate(self.t) + for i_wave in range(len(sub.waveforms)): + ind = np.logical_and((self.t-offset[i_wave] >= time_range[0]), (self.t-offset[i_wave] <= time_range[1])) + cc = signal.correlate(self.waveforms[i_wave][ind], sub.waveforms[i_wave][ind]) + dt = (np.argmax(cc) - (len(self.waveforms[i_wave][ind]) - 1)) * (self.t[1] - self.t[0]) + shift_list.append(dt) + return shift_list + + def interpolate(self, t_new): + """ + not in-place + """ + sec_new = WaveformSection() + sec_new.coord_list = self.coord_list + sec_new.t = np.copy(t_new) + sec_new.waveforms = [] + for wave in self.waveforms: + f = interp1d(self.t, wave, fill_value="extrapolate") + sec_new.waveforms.append(f(t_new)) + return sec_new + + def filter(self, freq_range, verbose=1): + if (not isinstance(freq_range, tuple)): + sys.exit(f"{freq_range} is not a tuple\n") + if (len(freq_range)!=2): + sys.exit(f"incorrect frequency range {freq_range}\n") + if (verbose==1): print(f'apply filter {freq_range}\n') + s = self.t[1] - self.t[0] + if (not (freq_range[0] >= 0.0)): + sos = signal.butter(4, freq_range[1] * s * 2, 'lowpass', output='sos') + else: + sos = signal.butter(4, [freq_range[0] * s * 2, freq_range[1] * s * 2], 'bandpass', output='sos') + for i_wave in range(len(self.waveforms)): + wave = self.waveforms[i_wave] + self.waveforms[i_wave] = signal.sosfiltfilt(sos, wave, padtype=None) + + def get_max_time(self, time_range=None): + if time_range is None: time_range = (self.t[0], self.t[-1]) + max_time_list = [] + max_amp_list = [] + for i_wave in range(len(self.waveforms)): + wave = self.waveforms[i_wave] + ind = np.logical_and((self.t >= time_range[0]), (self.t <= time_range[1])) + wave_abs = np.absolute(wave) + ind_max = np.argmax(wave_abs[ind]) + max_amp_list.append(np.amax(wave_abs[ind])) + max_time_list.append(time_range[0] + ind_max * (self.t[1] - self.t[0])) + return max_time_list, max_amp_list + + def get_first_arrival_time(self, threshold=0.1, time_range=None): + if time_range is None: time_range = (self.t[0], self.t[-1]) + max_time_list, max_amp_list = self.get_max_time() + first_arrival_time_list = [] + for i_wave in range(len(self.waveforms)): + wave = self.waveforms[i_wave] + ind = np.logical_and((self.t >= time_range[0]), (self.t <= time_range[1])) + wave = wave[ind] + wave_abs = np.absolute(wave) + for it in range(1, len(wave)-1): + if ((wave_abs[it] >= wave_abs[it-1]) and (wave_abs[it] >= wave_abs[it+1]) and (wave_abs[it]>=threshold*max_amp_list[i_wave])): + arrival = time_range[0] + it * (self.t[1] - self.t[0]) + break + first_arrival_time_list.append(arrival) + return first_arrival_time_list + + + def plot_waveforms_fill(self, *args, fig_num=0, time_range=None, offset=None, + x_axis_trans=(1.0, 0.0, 0.0, 0.0), + fill=True, fill_color=('blue', 'red'), + norm_with=None, normalize=1.0, label=None, **kwargs): + fig = plt.figure(num=fig_num) # pull out the existed figure + if offset is None: offset = [0.0] * len(self.waveforms) + if time_range is None: time_range = (self.t[0], self.t[-1]) + if norm_with is None: norm_with = self + max_val = 0.0 + for i_wave in range(len(self.waveforms)): + wave = norm_with.waveforms[i_wave] + ind = np.logical_and((norm_with.t-offset[i_wave] >= time_range[0]), (norm_with.t-offset[i_wave] <= time_range[1])) + if (np.amax(np.absolute(wave)) > max_val): + max_val = np.amax(np.absolute(wave[ind])) + for i_wave in range(len(self.waveforms)): + x_axis_val = x_axis_trans[0] * self.coord_list[i_wave][0] + \ + x_axis_trans[1] * self.coord_list[i_wave][1] + \ + x_axis_trans[2] * self.coord_list[i_wave][2] + \ + x_axis_trans[3] + ind = np.logical_and((self.t-offset[i_wave] >= time_range[0]), (self.t-offset[i_wave] <= time_range[1])) + wave_norm = self.waveforms[i_wave][ind] / max_val * normalize + x_axis_val + t_ind = self.t[ind] - offset[i_wave] + if ((i_wave == 0) and (label is not None)): + plt.plot(wave_norm, t_ind, *args, label=label, **kwargs) + else: + plt.plot(wave_norm, t_ind, *args, **kwargs) + if fill: + ax = plt.gca() + zero = np.zeros_like(wave_norm) + x_axis_val + ax.fill_betweenx(t_ind, wave_norm, zero, where=wave_norm >= zero , facecolor=fill_color[0]) + ax.fill_betweenx(t_ind, wave_norm, zero, where=wave_norm <= zero , facecolor=fill_color[1]) + + def plot_waveforms_ax_fill(self, *args, ax=None, time_range=None, offset=None, + x_axis_trans=(1.0, 0.0, 0.0, 0.0), + is_time_axis_x = False, + fill=True, fill_color=('blue', 'red'), + norm_with=None, normalize=1.0, label=None, **kwargs): + #fig = plt.figure(num=fig_num) # pull out the existed figure + if ax is None: ax = plt.gca() + if offset is None: offset = [0.0] * len(self.waveforms) + if time_range is None: time_range = (self.t[0], self.t[-1]) + if norm_with is None: norm_with = self + max_val = 0.0 + for i_wave in range(len(self.waveforms)): + wave = norm_with.waveforms[i_wave] + ind = np.logical_and((norm_with.t-offset[i_wave] >= time_range[0]), (norm_with.t-offset[i_wave] <= time_range[1])) + if (np.amax(np.absolute(wave)) > max_val): + max_val = np.amax(np.absolute(wave[ind])) + for i_wave in range(len(self.waveforms)): + x_axis_val = x_axis_trans[0] * self.coord_list[i_wave][0] + \ + x_axis_trans[1] * self.coord_list[i_wave][1] + \ + x_axis_trans[2] * self.coord_list[i_wave][2] + \ + x_axis_trans[3] + ind = np.logical_and((self.t-offset[i_wave] >= time_range[0]), (self.t-offset[i_wave] <= time_range[1])) + wave_norm = self.waveforms[i_wave][ind] / max_val * normalize + x_axis_val + t_ind = self.t[ind] - offset[i_wave] + if ((i_wave == 0) and (label is not None)): + if is_time_axis_x: + ax.plot(t_ind, wave_norm, *args, label=label, **kwargs) + else: + ax.plot(wave_norm, t_ind, *args, label=label, **kwargs) + else: + if is_time_axis_x: + ax.plot(t_ind, wave_norm, *args, **kwargs) + else: + ax.plot(wave_norm, t_ind, *args, **kwargs) + if fill: + #ax = plt.gca() + zero = np.zeros_like(wave_norm) + x_axis_val + if is_time_axis_x: + ax.fill_between(t_ind, wave_norm, zero, where=wave_norm >= zero , facecolor=fill_color[0]) + ax.fill_between(t_ind, wave_norm, zero, where=wave_norm <= zero , facecolor=fill_color[1]) + else: + ax.fill_betweenx(t_ind, wave_norm, zero, where=wave_norm >= zero , facecolor=fill_color[0]) + ax.fill_betweenx(t_ind, wave_norm, zero, where=wave_norm <= zero , facecolor=fill_color[1]) + + def get_waveforms_from_specfem(self, fn_list, fn_stations, + t=None, tstart=None, dt=None, nt=None, + verbose=1): + if (verbose==1): print('plotting waveforms from SPECFEM results\n') + n_stations = len(fn_list) + if (verbose==1): print(f'there are {n_stations} waveforms\n') + if (verbose==1): print(f'reading station information from {fn_stations}\n') + try: + coord_list = [None] * n_stations + with open(fn_stations, 'r') as f_stations: + lines = f_stations.readlines() + for line in lines: + line_segs = split_line(line) + nt_name = line_segs[0] + sta_name = line_segs[1] + sta_lat = float(line_segs[2]) + sta_lon = float(line_segs[3]) + sta_depth = float(line_segs[5]) + for i_fn in range(len(fn_list)): + fn = os.path.split(fn_list[i_fn])[-1] + if (fn.startswith(f'{nt_name}.{sta_name}') or fn.startswith(f'{sta_name}.{nt_name}')): + if coord_list[i_fn] is not None: + raise(ValueError('file list conflict')) + coord_list[i_fn] = (sta_lon / 1000.0, sta_lat / 1000.0, sta_depth / 1000.0) + break + for i_fn in range(len(fn_list)): + if coord_list[i_fn] is None: + raise(ValueError(f'station does not exist for {fn_list[i_fn]}')) + except Exception as e: + sys.exit(f"{e}\n") + self.coord_list = coord_list + + waveforms = [] + if t is not None: + if (verbose==1): print(f'use time t0={t[0]}, dt={t[1]-t[0]}, nt={nt}\n') + elif (tstart is not None) and (dt is not None) and (nt is not None): + try: + t = np.arange(0, nt, dtype=float) *dt + tstart + except Exception as e: + sys.exit(f"{e}\n") + if (verbose==1): print(f'use time t0={t[0]}, dt={t[1]-t[0]}, nt={nt}\n') + else: + if (verbose==1): print('find time in waveform files\n') + for i_fn in range(len(fn_list)): + try: + if (verbose==1): print(f'reading waveform file {fn_list[i_fn]}\n') + wave = np.loadtxt(fn_list[i_fn]) + if t is None: + if (tstart is None): tstart = wave[0,0] + if (dt is None): dt = wave[1,0] - wave[0,0] + if (nt is None): nt = len(wave) + try: + t = np.arange(0, nt, dtype=float) *dt + tstart + except Exception as e: + sys.exit(f"{e}\n") + if (verbose==1): print(f'use time t0={t[0]}, dt={t[1]-t[0]}, nt={nt}\n') + if (wave.shape[0] != len(t)): + raise (ValueError(f'file length inconsistant: {fn_list[i_fn]}')) + waveforms.append(wave[:,1]) + except Exception as e: + sys.exit(f"{e}\n") + self.waveforms = waveforms + self.t = t + + diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/CMTSOLUTION b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/CMTSOLUTION new file mode 100644 index 000000000..07356239f --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/CMTSOLUTION @@ -0,0 +1,13 @@ +PDE 1999 01 01 00 00 00.00 67000 67000 -25000 4.2 4.2 hom_explosion +event name: hom_explosion +time shift: 0.0000 +half duration: 5.0 +latorUTM: 33.6 +longorUTM: -118.4 +depth: 10.0 +Mrr: 1.000000e+23 +Mtt: 1.000000e+23 +Mpp: 1.000000e+23 +Mrt: 0.000000 +Mrp: 0.000000 +Mtp: 0.000000 diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/FORCESOLUTION b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/FORCESOLUTION new file mode 100644 index 000000000..e69de29bb diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/Par_file b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/Par_file new file mode 100644 index 000000000..c8d6b9b01 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/Par_file @@ -0,0 +1,405 @@ +#----------------------------------------------------------- +# +# Simulation input parameters +# +#----------------------------------------------------------- + +# forward or adjoint simulation +# 1 = forward, 2 = adjoint, 3 = both simultaneously +SIMULATION_TYPE = 1 +# 0 = earthquake simulation, 1/2/3 = three steps in noise simulation +NOISE_TOMOGRAPHY = 0 +SAVE_FORWARD = .false. + +# solve a full FWI inverse problem from a single calling program with no I/Os, storing everything in memory, +# or run a classical forward or adjoint problem only and save the seismograms and/or sensitivity kernels to disk (with costlier I/Os) +INVERSE_FWI_FULL_PROBLEM = .false. + +# UTM projection parameters +# Use a negative zone number for the Southern hemisphere: +# The Northern hemisphere corresponds to zones +1 to +60, +# The Southern hemisphere corresponds to zones -1 to -60. +UTM_PROJECTION_ZONE = 11 +SUPPRESS_UTM_PROJECTION = .true. + +# number of MPI processors +NPROC = 4 + +# time step parameters +NSTEP = 3000 +DT = 0.02 + +# set to true to use local-time stepping (LTS) +LTS_MODE = .false. + +# Partitioning algorithm for decompose_mesh +# choose partitioner: 1==SCOTCH (default), 2==METIS, 3==PATOH, 4==ROWS_PART +PARTITIONING_TYPE = 1 + +#----------------------------------------------------------- +# +# LDDRK time scheme +# +#----------------------------------------------------------- +USE_LDDRK = .false. +INCREASE_CFL_FOR_LDDRK = .false. +RATIO_BY_WHICH_TO_INCREASE_IT = 1.4 + +#----------------------------------------------------------- +# +# Mesh +# +#----------------------------------------------------------- + +# Number of nodes for 2D and 3D shape functions for hexahedra. +# We use either 8-node mesh elements (bricks) or 27-node elements. +# If you use our internal mesher, the only option is 8-node bricks (27-node elements are not supported). +NGNOD = 8 + +# models: +# available options are: +# default (model parameters described by mesh properties) +# 1D models available are: +# 1d_prem,1d_socal,1d_cascadia +# 3D models available are: +# aniso,external,gll,salton_trough,tomo,SEP,coupled,... +MODEL = default + +# path for external tomographic models files +TOMOGRAPHY_PATH = DATA/tomo_files/ +# if you are using a SEP model (oil-industry format) +SEP_MODEL_DIRECTORY = DATA/my_SEP_model/ + +#----------------------------------------------------------- + +# parameters describing the model +APPROXIMATE_OCEAN_LOAD = .false. +TOPOGRAPHY = .false. +ATTENUATION = .false. +ANISOTROPY = .false. +GRAVITY = .false. + +# in case of attenuation, reference frequency in Hz at which the velocity values in the velocity model are given (unused otherwise) +ATTENUATION_f0_REFERENCE = 18.d0 + +# attenuation period range over which we try to mimic a constant Q factor +MIN_ATTENUATION_PERIOD = 999999998.d0 +MAX_ATTENUATION_PERIOD = 999999999.d0 +# ignore this range and ask the code to compute it automatically instead based on the estimated resolution of the mesh (use this unless you know what you are doing) +COMPUTE_FREQ_BAND_AUTOMATIC = .true. + +# Olsen's constant for Q_mu = constant * V_s attenuation rule +USE_OLSEN_ATTENUATION = .false. +OLSEN_ATTENUATION_RATIO = 0.05 + +#----------------------------------------------------------- +# +# Absorbing boundary conditions +# +#----------------------------------------------------------- + +# C-PML boundary conditions for a regional simulation +# (if set to .false., and STACEY_ABSORBING_CONDITIONS is also set to .false., you get a free surface instead +# in the case of elastic or viscoelastic mesh elements, and a rigid surface in the case of acoustic (fluid) elements +PML_CONDITIONS = .true. + +# C-PML top surface +PML_INSTEAD_OF_FREE_SURFACE = .false. + +# C-PML dominant frequency +f0_FOR_PML = 1.0 + +# parameters used to rotate C-PML boundary conditions by a given angle (not completed yet) +# ROTATE_PML_ACTIVATE = .false. +# ROTATE_PML_ANGLE = 0. + +# absorbing boundary conditions for a regional simulation +# (if set to .false., and PML_CONDITIONS is also set to .false., you get a free surface instead +# in the case of elastic or viscoelastic mesh elements, and a rigid surface in the case of acoustic (fluid) elements +STACEY_ABSORBING_CONDITIONS = .false. + +# absorbing top surface (defined in mesh as 'free_surface_file') +STACEY_INSTEAD_OF_FREE_SURFACE = .false. + +# When STACEY_ABSORBING_CONDITIONS is set to .true. : +# absorbing conditions are defined in xmin, xmax, ymin, ymax and zmin +# this option BOTTOM_FREE_SURFACE can be set to .true. to +# make zmin free surface instead of absorbing condition +BOTTOM_FREE_SURFACE = .false. + +#----------------------------------------------------------- +# +# undoing attenuation and/or PMLs for sensitivity kernel calculations +# +#----------------------------------------------------------- + +# to undo attenuation and/or PMLs for sensitivity kernel calculations or forward runs with SAVE_FORWARD +# use the flag below. It performs undoing of attenuation and/or of PMLs in an exact way for sensitivity kernel calculations +# but requires disk space for temporary storage, and uses a significant amount of memory used as buffers for temporary storage. +# When that option is on the second parameter indicates how often the code dumps restart files to disk (if in doubt, use something between 100 and 1000). +UNDO_ATTENUATION_AND_OR_PML = .false. +NT_DUMP_ATTENUATION = 500 + +#----------------------------------------------------------- +# +# Visualization +# +#----------------------------------------------------------- + +# save AVS or OpenDX movies +# MOVIE_TYPE = 1 to show the top surface +# MOVIE_TYPE = 2 to show all the external faces of the mesh +CREATE_SHAKEMAP = .false. +MOVIE_SURFACE = .false. +MOVIE_TYPE = 1 +MOVIE_VOLUME = .false. +SAVE_DISPLACEMENT = .false. +MOVIE_VOLUME_STRESS = .false. +USE_HIGHRES_FOR_MOVIES = .false. +NTSTEP_BETWEEN_FRAMES = 200 +HDUR_MOVIE = 0.0 + +# save AVS or OpenDX mesh files to check the mesh +SAVE_MESH_FILES = .true. + +# path to store the local database file on each node +LOCAL_PATH = ./DATABASES_MPI + +# interval at which we output time step info and max of norm of displacement +NTSTEP_BETWEEN_OUTPUT_INFO = 500 + +#----------------------------------------------------------- +# +# Sources +# +#----------------------------------------------------------- + +# sources and receivers Z coordinates given directly (i.e. as their true position) instead of as their depth +USE_SOURCES_RECEIVERS_Z = .true. + +# use a (tilted) FORCESOLUTION force point source (or several) instead of a CMTSOLUTION moment-tensor source. +# This can be useful e.g. for oil industry foothills simulations or asteroid simulations +# in which the source is a vertical force, normal force, tilted force, impact etc. +# If this flag is turned on, the FORCESOLUTION file must be edited by giving: +# - the corresponding time-shift parameter, +# - the half duration (hdur, in s) for Gaussian/Step function, dominant frequency (f0, in Hz) for Ricker, +# - the coordinates of the source, +# - the source time function type (0=Gaussian function, 1=Ricker wavelet, 2=Step function), +# - the magnitude of the force source, +# - the components of a (non necessarily unitary) direction vector for the force source in the E/N/Z_UP basis. +# The direction vector is made unitary internally in the code and thus only its direction matters here; +# its norm is ignored and the norm of the force used is the factor force source times the source time function. +USE_FORCE_POINT_SOURCE = .true. + +# set to true to use a Ricker source time function instead of the source time functions set by default +# to represent a (tilted) FORCESOLUTION force point source or a CMTSOLUTION moment-tensor source. +USE_RICKER_TIME_FUNCTION = .true. + +# use an external source time function +# you must add a file with your source time function and the file name path +# relative to working directory at the end of CMTSOLUTION or FORCESOLUTION file +# (with multiple sources, one file per source is required). +# This file must have a single column containing the amplitudes of the source at all time steps; +# time step size used must be equal to DT as defined at the beginning of this Par_file. +# Be sure when this option is .false. to remove the name of stf file in CMTSOLUTION or FORCESOLUTION +USE_EXTERNAL_SOURCE_FILE = .false. + +# print source time function +PRINT_SOURCE_TIME_FUNCTION = .true. + +# source encoding +# (for acoustic simulations only for now) determines source encoding factor +/-1 depending on sign of moment tensor +# (see e.g. Krebs et al., 2009. Fast full-wavefield seismic inversion using encoded sources, Geophysics, 74 (6), WCC177-WCC188.) +USE_SOURCE_ENCODING = .false. + +#----------------------------------------------------------- +# +# Seismograms +# +#----------------------------------------------------------- + +# interval in time steps for writing of seismograms +NTSTEP_BETWEEN_OUTPUT_SEISMOS = 10000 + +# set to n to reduce the sampling rate of output seismograms by a factor of n +# defaults to 1, which means no down-sampling +NTSTEP_BETWEEN_OUTPUT_SAMPLE = 1 + +# decide if we save displacement, velocity, acceleration and/or pressure in forward runs (they can be set to true simultaneously) +# currently pressure seismograms are implemented in acoustic (i.e. fluid) elements only +SAVE_SEISMOGRAMS_DISPLACEMENT = .true. +SAVE_SEISMOGRAMS_VELOCITY = .false. +SAVE_SEISMOGRAMS_ACCELERATION = .false. +SAVE_SEISMOGRAMS_PRESSURE = .false. # currently implemented in acoustic (i.e. fluid) elements only + +# option to save strain seismograms +# this option is useful for strain Green's tensor +SAVE_SEISMOGRAMS_STRAIN = .false. + +# save seismograms also when running the adjoint runs for an inverse problem +# (usually they are unused and not very meaningful, leave this off in almost all cases) +SAVE_SEISMOGRAMS_IN_ADJOINT_RUN = .false. + +# save seismograms in binary or ASCII format (binary is smaller but may not be portable between machines) +USE_BINARY_FOR_SEISMOGRAMS = .false. + +# output seismograms in Seismic Unix format (binary with 240-byte-headers) +SU_FORMAT = .false. + +# output seismograms in ASDF (requires asdf-library) +ASDF_FORMAT = .false. + +# output seismograms in HDF5 (requires hdf5-library and WRITE_SEISMOGRAMS_BY_MAIN) +HDF5_FORMAT = .false. + +# decide if main process writes all the seismograms or if all processes do it in parallel +WRITE_SEISMOGRAMS_BY_MAIN = .false. + +# save all seismograms in one large combined file instead of one file per seismogram +# to avoid overloading shared non-local file systems such as LUSTRE or GPFS for instance +SAVE_ALL_SEISMOS_IN_ONE_FILE = .false. + +# use a trick to increase accuracy of pressure seismograms in fluid (acoustic) elements: +# use the second derivative of the source for the source time function instead of the source itself, +# and then record -potential_acoustic() as pressure seismograms instead of -potential_dot_dot_acoustic(); +# this is mathematically equivalent, but numerically significantly more accurate because in the explicit +# Newmark time scheme acceleration is accurate at zeroth order while displacement is accurate at second order, +# thus in fluid elements potential_dot_dot_acoustic() is accurate at zeroth order while potential_acoustic() +# is accurate at second order and thus contains significantly less numerical noise. +USE_TRICK_FOR_BETTER_PRESSURE = .false. + +#----------------------------------------------------------- +# +# Energy calculation +# +#----------------------------------------------------------- + +# to plot energy curves, for instance to monitor how CPML absorbing layers behave; +# should be turned OFF in most cases because a bit expensive +OUTPUT_ENERGY = .false. +# every how many time steps we compute energy (which is a bit expensive to compute) +NTSTEP_BETWEEN_OUTPUT_ENERGY = 10 + +#----------------------------------------------------------- +# +# Adjoint kernel outputs +# +#----------------------------------------------------------- + +# interval in time steps for reading adjoint traces +# 0 = read the whole adjoint sources at start time +NTSTEP_BETWEEN_READ_ADJSRC = 0 + +# read adjoint sources using ASDF (requires asdf-library) +READ_ADJSRC_ASDF = .false. + +# this parameter must be set to .true. to compute anisotropic kernels +# in crust and mantle (related to the 21 Cij in geographical coordinates) +# default is .false. to compute isotropic kernels (related to alpha and beta) +ANISOTROPIC_KL = .false. + +# compute transverse isotropic kernels (alpha_v,alpha_h,beta_v,beta_h,eta,rho) +# rather than fully anisotropic kernels in case ANISOTROPIC_KL is set to .true. +SAVE_TRANSVERSE_KL = .false. + +# this parameter must be set to .true. to compute anisotropic kernels for +# cost function using velocity observable rather than displacement +ANISOTROPIC_VELOCITY_KL = .false. + +# outputs approximate Hessian for preconditioning +APPROXIMATE_HESS_KL = .false. + +# save Moho mesh and compute Moho boundary kernels +SAVE_MOHO_MESH = .false. + +#----------------------------------------------------------- +# +# Coupling with an injection technique (DSM, AxiSEM, or FK) +# +#----------------------------------------------------------- +COUPLE_WITH_INJECTION_TECHNIQUE = .false. +INJECTION_TECHNIQUE_TYPE = 3 # 1 = DSM, 2 = AxiSEM, 3 = FK +MESH_A_CHUNK_OF_THE_EARTH = .false. +TRACTION_PATH = DATA/AxiSEM_tractions/3/ +FKMODEL_FILE = FKmodel +RECIPROCITY_AND_KH_INTEGRAL = .false. # does not work yet + +#----------------------------------------------------------- +## +## Prescribed wavefield discontinuity on an interface +## +##----------------------------------------------------------- +IS_WAVEFIELD_DISCONTINUITY = .true. + +#----------------------------------------------------------- +# +# Run modes +# +#----------------------------------------------------------- + +# Simultaneous runs +# added the ability to run several calculations (several earthquakes) +# in an embarrassingly-parallel fashion from within the same run; +# this can be useful when using a very large supercomputer to compute +# many earthquakes in a catalog, in which case it can be better from +# a batch job submission point of view to start fewer and much larger jobs, +# each of them computing several earthquakes in parallel. +# To turn that option on, set parameter NUMBER_OF_SIMULTANEOUS_RUNS to a value greater than 1. +# To implement that, we create NUMBER_OF_SIMULTANEOUS_RUNS MPI sub-communicators, +# each of them being labeled "my_local_mpi_comm_world", and we use them +# in all the routines in "src/shared/parallel.f90", except in MPI_ABORT() because in that case +# we need to kill the entire run. +# When that option is on, of course the number of processor cores used to start +# the code in the batch system must be a multiple of NUMBER_OF_SIMULTANEOUS_RUNS, +# all the individual runs must use the same number of processor cores, +# which as usual is NPROC in the Par_file, +# and thus the total number of processor cores to request from the batch system +# should be NUMBER_OF_SIMULTANEOUS_RUNS * NPROC. +# All the runs to perform must be placed in directories called run0001, run0002, run0003 and so on +# (with exactly four digits). +# +# Imagine you have 10 independent calculations to do, each of them on 100 cores; you have three options: +# +# 1/ submit 10 jobs to the batch system +# +# 2/ submit a single job on 1000 cores to the batch, and in that script create a sub-array of jobs to start 10 jobs, +# each running on 100 cores (see e.g. http://www.schedmd.com/slurmdocs/job_array.html ) +# +# 3/ submit a single job on 1000 cores to the batch, start SPECFEM3D on 1000 cores, create 10 sub-communicators, +# cd into one of 10 subdirectories (called e.g. run0001, run0002,... run0010) depending on the sub-communicator +# your MPI rank belongs to, and run normally on 100 cores using that sub-communicator. +# +# The option below implements 3/. +# +NUMBER_OF_SIMULTANEOUS_RUNS = 1 + +# if we perform simultaneous runs in parallel, if only the source and receivers vary between these runs +# but not the mesh nor the model (velocity and density) then we can also read the mesh and model files +# from a single run in the beginning and broadcast them to all the others; for a large number of simultaneous +# runs for instance when solving inverse problems iteratively this can DRASTICALLY reduce I/Os to disk in the solver +# (by a factor equal to NUMBER_OF_SIMULTANEOUS_RUNS), and reducing I/Os is crucial in the case of huge runs. +# Thus, always set this option to .true. if the mesh and the model are the same for all simultaneous runs. +# In that case there is no need to duplicate the mesh and model file database (the content of the DATABASES_MPI +# directories) in each of the run0001, run0002,... directories, it is sufficient to have one in run0001 +# and the code will broadcast it to the others) +BROADCAST_SAME_MESH_AND_MODEL = .true. + +#----------------------------------------------------------- + +# set to true to use GPUs +GPU_MODE = .false. + +# ADIOS Options for I/Os +ADIOS_ENABLED = .false. +ADIOS_FOR_DATABASES = .false. +ADIOS_FOR_MESH = .false. +ADIOS_FOR_FORWARD_ARRAYS = .false. +ADIOS_FOR_KERNELS = .false. +ADIOS_FOR_UNDO_ATTENUATION = .false. + +# HDF5 Database I/O +# (note the flags for HDF5 and ADIOS are mutually exclusive, only one can be used) +HDF5_ENABLED = .false. +HDF5_FOR_MOVIES = .false. +HDF5_IO_NODES = 0 # HDF5 IO server with number of IO dedicated procs + diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/STATIONS b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/STATIONS new file mode 100644 index 000000000..f9175a451 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/STATIONS @@ -0,0 +1,4 @@ +TS01 TS 25000.0 10000.0 0.0 0.0 +TS02 TS 25000.0 20000.0 0.0 0.0 +TS03 TS 25000.0 30000.0 0.0 0.0 +TS04 TS 25000.0 40000.0 0.0 0.0 diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/STATIONS_FILTERED b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/STATIONS_FILTERED new file mode 100644 index 000000000..6b1e6c492 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/STATIONS_FILTERED @@ -0,0 +1,4 @@ + TS01 TS 25000.000000000000 10000.000000000000 0.000000000000 0.000000000000 + TS02 TS 25000.000000000000 20000.000000000000 0.000000000000 0.000000000000 + TS03 TS 25000.000000000000 30000.000000000000 0.000000000000 0.000000000000 + TS04 TS 25000.000000000000 40000.000000000000 0.000000000000 0.000000000000 diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/meshfem3D_files/Mesh_Par_file b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/meshfem3D_files/Mesh_Par_file new file mode 100644 index 000000000..2610f25ff --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/meshfem3D_files/Mesh_Par_file @@ -0,0 +1,103 @@ +#----------------------------------------------------------- +# +# Meshing input parameters +# +#----------------------------------------------------------- + +# coordinates of mesh block in latitude/longitude and depth in km +LATITUDE_MIN = -10000.0 +LATITUDE_MAX = 60000.0 +LONGITUDE_MIN = -10000.0 +LONGITUDE_MAX = 60000.0 +DEPTH_BLOCK_KM = 60.0 +UTM_PROJECTION_ZONE = 11 +SUPPRESS_UTM_PROJECTION = .true. + +# file that contains the interfaces of the model / mesh +INTERFACES_FILE = interfaces.dat + +# file that contains the cavity +CAVITY_FILE = no_cavity.dat + +# number of elements at the surface along edges of the mesh at the surface +# (must be 8 * multiple of NPROC below if mesh is not regular and contains mesh doublings) +# (must be multiple of NPROC below if mesh is regular) +NEX_XI = 28 +NEX_ETA = 28 + +# number of MPI processors along xi and eta (can be different) +NPROC_XI = 1 +NPROC_ETA = 1 + +#----------------------------------------------------------- +# +# Doubling layers +# +#----------------------------------------------------------- + +# Regular/irregular mesh +USE_REGULAR_MESH = .true. +# Only for irregular meshes, number of doubling layers and their position +NDOUBLINGS = 0 +# NZ_DOUBLING_1 is the parameter to set up if there is only one doubling layer +# (more doubling entries can be added if needed to match NDOUBLINGS value) +NZ_DOUBLING_1 = 11 +NZ_DOUBLING_2 = 0 + +#----------------------------------------------------------- +# +# Visualization +# +#----------------------------------------------------------- + +# create mesh files for visualisation or further checking +CREATE_ABAQUS_FILES = .false. +CREATE_DX_FILES = .false. +CREATE_VTK_FILES = .false. + +# stores mesh files as cubit-exported files into directory MESH/ (for single process run) +SAVE_MESH_AS_CUBIT = .true. + +# path to store the databases files +LOCAL_PATH = ./DATABASES_MPI + +#----------------------------------------------------------- +# +# CPML +# +#----------------------------------------------------------- + +# CPML perfectly matched absorbing layers +THICKNESS_OF_X_PML = 10000.0 +THICKNESS_OF_Y_PML = 10000.0 +THICKNESS_OF_Z_PML = 10000.0 + +#----------------------------------------------------------- +# +# Domain materials +# +#----------------------------------------------------------- + +# number of materials +NMATERIALS = 2 +# define the different materials in the model as: +# #material_id #rho #vp #vs #Q_Kappa #Q_mu #anisotropy_flag #domain_id +# Q_Kappa : Q_Kappa attenuation quality factor +# Q_mu : Q_mu attenuation quality factor +# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90 +# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic +1 3380.0 8080.0 4485.0 9999. 40.0 0 2 +2 2600.0 5800.0 3198.0 9999. 40.0 0 2 + +#----------------------------------------------------------- +# +# Domain regions +# +#----------------------------------------------------------- + +# number of regions +NREGIONS = 2 +# define the different regions of the model as : +#NEX_XI_BEGIN #NEX_XI_END #NEX_ETA_BEGIN #NEX_ETA_END #NZ_BEGIN #NZ_END #material_id +1 28 1 28 1 10 1 +1 28 1 28 11 24 2 diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/meshfem3D_files/interfaces.dat b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/meshfem3D_files/interfaces.dat new file mode 100644 index 000000000..2d1101962 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/meshfem3D_files/interfaces.dat @@ -0,0 +1,18 @@ +# number of interfaces + 2 +# +# We describe each interface below, structured as a 2D-grid, with several parameters : +# number of points along XI and ETA, minimal XI ETA coordinates +# and spacing between points which must be constant. +# Then the records contain the Z coordinates of the NXI x NETA points. +# +# interface number 1 +# SUPPRESS_UTM_PROJECTION NXI NETA LONG_MIN LAT_MIN SPACING_XI SPACING_ETA + .true. 2 2 -10000.0 -10000.0 70000.0 70000.0 + moho.dat + .true. 2 2 -10000.0 -10000.0 70000.0 70000.0 + top.dat +# for each layer, we give the number of spectral elements in the vertical direction +# layer number 1 (bottom layer) + 10 + 14 diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/meshfem3D_files/moho.dat b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/meshfem3D_files/moho.dat new file mode 100644 index 000000000..01cb63a8c --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/meshfem3D_files/moho.dat @@ -0,0 +1,4 @@ +-35000.0 +-35000.0 +-35000.0 +-35000.0 diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/meshfem3D_files/top.dat b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/meshfem3D_files/top.dat new file mode 100644 index 000000000..68ea011c5 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/DATA/meshfem3D_files/top.dat @@ -0,0 +1,4 @@ +0.0 +0.0 +0.0 +0.0 diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/fk_model b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/fk_model new file mode 100644 index 000000000..7520f3074 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/fk_model @@ -0,0 +1,14 @@ +NLAYER 2 +LAYER 1 2600.0 5800.0 3198.0 0.0 +LAYER 2 3380.0 8080.0 4485.0 -35000.0 +INCIDENT_WAVE P +BACK_AZIMUTH 270.0 +TAKE_OFF 15.0 +ORIGIN_WAVEFRONT 0.0 25000.0 -75000.0 +ORIGIN_TIME 0.0 +NSTEP 3000 +deltat 0.02 +FREQUENCY_MAX 1.0 +FREQUENCY_SAMPLING 32.0 +TIME_WINDOW 400.0 +AMPLITUDE 1.0 diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/plot_seismograms.py b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/plot_seismograms.py new file mode 100644 index 000000000..cd565d847 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/plot_seismograms.py @@ -0,0 +1,53 @@ +import sys +sys.path.insert(1, '../plot') +from plot_utils_injection import WaveformSection +import matplotlib.pyplot as plt +import os + +Tmin = 1.5 +Tmax = 50.0 + +fig = plt.figure(figsize=(12, 9)) + +wave_dir = 'OUTPUT_FILES' + +sta_list = [f'TS{i:02d}' for i in range(1,5)] + +fn_list = [os.path.join(wave_dir, 'TS.'+sta+'.Z.fkd') for sta in sta_list] + +sec_fk = WaveformSection() + +sec_fk.get_waveforms_from_specfem(fn_list=fn_list, fn_stations=os.path.join(wave_dir, 'STATIONS')) + + +wave_dir = 'OUTPUT_FILES' + +sta_list = [f'TS{i:02d}' for i in range(1,5)] + +fn_list = [os.path.join(wave_dir, 'TS.'+sta+'.BXZ.semd') for sta in sta_list] + +sec = WaveformSection() + +sec.get_waveforms_from_specfem(fn_list=fn_list, fn_stations=os.path.join(wave_dir, 'STATIONS')) + +sec_fk.filter((-1.0, 1.0/Tmin)) +sec.filter((-1.0, 1.0/Tmin)) + +sec_diff = sec - sec_fk + +offset = sec_fk.get_first_arrival_time() +print(offset) + +print(sec_fk.difference_ratio(sec, offset=offset, time_range=(-10.0, 40.0))) +sec.plot_waveforms_fill('k', fig_num=fig.number, normalize=10.0, offset=offset, time_range=(-10.0, 40.0), fill=False, label="injection", linewidth=2) +sec_diff.plot_waveforms_fill('r--', fig_num=fig.number, normalize=100.0, offset=offset, norm_with=sec_fk, time_range=(-10.0, 40.0), fill=False, label="(injection-FK)X10", linewidth=2) +plt.ylim([-10.0, 40.0]) +plt.xlim([0.0, 50.0]) +plt.legend(prop={'size':20}, loc=4) +plt.xlabel('x (km)', fontsize=20) +plt.ylabel('t (s)', fontsize=20) +plt.gca().invert_yaxis() +plt.gca().tick_params(labelsize=20) +#plt.text(-30.0, -10.0, '(a)', fontsize=20) +plt.show() +#plt.savefig("crust_Z.pdf") diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/run_this_example.sh b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/run_this_example.sh new file mode 100755 index 000000000..9415f7134 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/run_this_example.sh @@ -0,0 +1,108 @@ +#!/bin/bash + +## load all necessary modules here +module load intel/2020u4 intelmpi/2020u4 + +## set number of processors +NPROC=4 + +currentdir=`pwd` + +## set the path to SPECFEM here +specfem_dir=~/specfem3d + + +echo +echo " running a small example for wavefield injection" +echo + + +BASEMPIDIR=`grep ^LOCAL_PATH DATA/Par_file | cut -d = -f 2 ` +mkdir -p OUTPUT_FILES/ +rm -rf OUTPUT_FILES/* +mkdir -p $BASEMPIDIR +rm -rf $BASEMPIDIR/* + +# links executables +mkdir -p bin +cd bin/ +rm -f * +ln -s ${specfem_dir}/bin/xmeshfem3D +ln -s ${specfem_dir}/bin/xdecompose_mesh +ln -s ${specfem_dir}/bin/xgenerate_databases +#ln -s ${specfem_dir}/bin/xcombine_vol_data_vtk +ln -s ${specfem_dir}/bin/xspecfem3D +cd ../ + +# stores setup +cp DATA/meshfem3D_files/Mesh_Par_file OUTPUT_FILES/ +cp DATA/Par_file OUTPUT_FILES/ +cp DATA/FORCESOLUTION OUTPUT_FILES/ +cp DATA/STATIONS OUTPUT_FILES/ + +# prepare to run the mesher +# at this point the meshfem3D files must be well-prepared, +# the Par_file should be set properly to turn on IS_WAVEFIELD_DISCONTINUITY, +# an empty FORCESOLUTION file should be present to indicate +# there is no external source, +# and a wavefield_discontinuity_box file to designate the interface +mkdir -p MESH/ + +###################################################################### +## running meshfem3D, here I adopt a strategy that first run meshfem3D +## in seriel, output the mesh in CUBIT format, and then decompose +## the mesh with decompose_mesh +sed -i "/^NPROC/c\NPROC = 1" DATA/Par_file +sed -i "/^NPROC_XI/c\NPROC_XI = 1" DATA/meshfem3D_files/Mesh_Par_file +sed -i "/^NPROC_ETA/c\NPROC_ETA = 1" DATA/meshfem3D_files/Mesh_Par_file +sed -i "/^SAVE_MESH_AS_CUBIT/c\SAVE_MESH_AS_CUBIT = .true." DATA/meshfem3D_files/Mesh_Par_file + +echo +echo " running mesher in seriel..." +echo +./bin/xmeshfem3D + +mkdir -p MESH-default +cp -f MESH/* MESH-default + +sed -i "/^NPROC/c\NPROC = ${NPROC}" DATA/Par_file + +echo +echo " decomposing mesh..." +echo +./bin/xdecompose_mesh $NPROC ./MESH-default $BASEMPIDIR + +###################################################################### +## note that this can be replaced by directly running meshfem3D +## in parallel +#sed -i "/^NPROC/c\NPROC = ${NPROC}" DATA/Par_file +#sed -i "/^NPROC_XI/c\NPROC_XI = 2" DATA/meshfem3D_files/Mesh_Par_file +#sed -i "/^NPROC_ETA/c\NPROC_ETA = 2" DATA/meshfem3D_files/Mesh_Par_file +#sed -i "/^SAVE_MESH_AS_CUBIT/c\SAVE_MESH_AS_CUBIT = .false." DATA/meshfem3D_files/Mesh_Par_file +# +#echo +#echo " running mesher in parallel..." +#echo +#mpirun -np ${NPROC} ./bin/xmeshfem3D +######################################################################## + +echo +echo " generate databases..." +echo +mpirun -np $NPROC ./bin/xgenerate_databases + +# generate injection wavefield by producing +# DATABASES_MPI/proc*_wavefield_discontinuity.bin files +mpirun -np $NPROC ../fk_coupling/compute_fk_injection_field + +echo +echo " launch solver..." +echo +mpirun -np $NPROC ./bin/xspecfem3D + +# compute reference seismograms using FK +../fk_coupling/compute_fk_receiver + +# checks exit code +if [[ $? -ne 0 ]]; then exit 1; fi + diff --git a/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/wavefield_discontinuity_box b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/wavefield_discontinuity_box new file mode 100644 index 000000000..bae672e44 --- /dev/null +++ b/EXAMPLES/applications/wavefield_discontinuity/test_small_injection/wavefield_discontinuity_box @@ -0,0 +1,8 @@ +.false. # IS_TOP_WAVEFIELD_DISCONTINUITY +.false. # IS_EXTRAPOLATION_MODE +0.0 # x_min +50000.0 # x_max +0.0 # y_min +50000.0 # y_max +-50000.0 # z_min +0.0 # z_max diff --git a/setup/constants.h.in b/setup/constants.h.in index 09b34e89a..3bfed6bfe 100644 --- a/setup/constants.h.in +++ b/setup/constants.h.in @@ -287,6 +287,28 @@ integer, parameter :: IIN_veloc_dsm = 51, IIN_tract_dsm = 52 integer, parameter :: IIN_displ_axisem = 54 +!!----------------------------------------------------------- +!! +!! Prescribed wavefield discontinuity on an interface +!! +!!----------------------------------------------------------- + + +! I/O for wavefield discontinuity + integer, parameter :: IFILE_WAVEFIELD_DISCONTINUITY = 527 + +! file name for ASCII wavefield discontinuity interface file + character(len=*), parameter :: FNAME_WAVEFIELD_DISCONTINUITY_INTERFACE = 'wavefield_discontinuity_interface' + +! file name for ASCII wavefield discontinuity box file + character(len=*), parameter :: FNAME_WAVEFIELD_DISCONTINUITY_BOX = 'wavefield_discontinuity_box' + +! file name for binary wavefield discontinuity mesh file + character(len=*), parameter :: FNAME_WAVEFIELD_DISCONTINUITY_MESH = 'wavefield_discontinuity_mesh.bin' + +! file name for binary wavefield discontinuity solver database file + character(len=*), parameter :: FNAME_WAVEFIELD_DISCONTINUITY_DATABASE = 'wavefield_discontinuity_database.bin' + !!----------------------------------------------------------- !! !! mesh optimization diff --git a/src/decompose_mesh/decompose_mesh_par.F90 b/src/decompose_mesh/decompose_mesh_par.F90 index df5853233..9735094de 100644 --- a/src/decompose_mesh/decompose_mesh_par.F90 +++ b/src/decompose_mesh/decompose_mesh_par.F90 @@ -39,6 +39,8 @@ module decompose_mesh_par acoustic_elastic_poro_load,mesh2dual_ncommonnodes, & build_glob2loc_elmnts,build_glob2loc_nodes,build_interfaces,poro_elastic_repartitioning,moho_surface_repartitioning + !! setting up wavefield discontinuity interface + use part_decompose_mesh, only: write_wavefield_discontinuity_database implicit none ! note: the poroelastic repartitioning routine to parallelize the poroelastic-elastic interface @@ -120,5 +122,23 @@ module decompose_mesh_par integer,dimension(:),allocatable :: num_ispec_level integer :: num_p_level + ! wavefield discontinuity + !! boundary of wavefield discontinuity, read from database file + integer :: nb_wd + + !! boundary_to_ispec_wd(nb_wd) + !! the element the boundary belongs to, read from database file + !! each point on the boundary belongs to two sides of the boundary + !! here the element must be on the inner side of the boundary + integer, dimension(:), allocatable :: boundary_to_ispec_wd + + !! side_wd(nb_wd) + !! integers specifying which side the boundary is in the element + !! read from database file + !! side_wd = 1--8: only one vertex is on the boundary + !! side_wd = 9--20: only one edge is on the boundary + !! side_wd = 21--26: one face is on the boundary + integer, dimension(:), allocatable :: side_wd + end module decompose_mesh_par diff --git a/src/decompose_mesh/part_decompose_mesh.F90 b/src/decompose_mesh/part_decompose_mesh.F90 index e0c8a095b..e1468dffe 100644 --- a/src/decompose_mesh/part_decompose_mesh.F90 +++ b/src/decompose_mesh/part_decompose_mesh.F90 @@ -608,6 +608,62 @@ subroutine write_material_props_database(IIN_database,count_def_mat, & end subroutine write_material_props_database + !------------------------------------------------- + ! Write elements on the inner side of the wavefield + ! discontinuity interface + !------------------------------------------------- + subroutine write_wavefield_discontinuity_database(iproc, outputpath_name, & + nb_wd, boundary_to_ispec_wd, side_wd, & + nspec, glob2loc_elmnts, part) + use constants, only: FNAME_WAVEFIELD_DISCONTINUITY_MESH, & + IFILE_WAVEFIELD_DISCONTINUITY, MAX_STRING_LEN + implicit none + integer, intent(in) :: iproc + integer, intent(in) :: nspec + integer, intent(in) :: nb_wd + integer, dimension(:), pointer :: glob2loc_elmnts + integer, dimension(1:nspec) :: part + integer, dimension(1:nb_wd), intent(in) :: boundary_to_ispec_wd, side_wd + integer :: ispec, iside, ib + integer :: local_nb_wd + integer, dimension(:), allocatable :: local_boundary_to_ispec_wd, & + local_side_wd + character(len=MAX_STRING_LEN), intent(in) :: outputpath_name + character(len=MAX_STRING_LEN) :: prname + write(prname, "('proc',i6.6,'_')") iproc + open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & + file=trim(outputpath_name)//'/'//trim(prname)//& + trim(FNAME_WAVEFIELD_DISCONTINUITY_MESH), & + form='unformatted', action='write') + local_nb_wd = 0 + do ib = 1, nb_wd + ispec = boundary_to_ispec_wd(ib) + iside = side_wd(ib) + if (part(ispec) == iproc) then + local_nb_wd = local_nb_wd + 1 + endif + enddo + print *, 'partition', iproc, ' has ', local_nb_wd, ' elements on wavefield discontinuity interface' + allocate(local_boundary_to_ispec_wd(local_nb_wd)) + allocate(local_side_wd(local_nb_wd)) + local_nb_wd = 0 + do ib = 1, nb_wd + ispec = boundary_to_ispec_wd(ib) + iside = side_wd(ib) + if (part(ispec) == iproc) then + local_nb_wd = local_nb_wd + 1 + local_boundary_to_ispec_wd(local_nb_wd) = glob2loc_elmnts(ispec-1)+1 + local_side_wd(local_nb_wd) = iside + endif + enddo + write(IFILE_WAVEFIELD_DISCONTINUITY) local_nb_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) local_boundary_to_ispec_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) local_side_wd + close(IFILE_WAVEFIELD_DISCONTINUITY) + deallocate(local_boundary_to_ispec_wd, local_side_wd) + end subroutine write_wavefield_discontinuity_database + + !-------------------------------------------------- ! Write elements on boundaries (and their four nodes on boundaries) ! pertaining to iproc partition in the corresponding Database diff --git a/src/decompose_mesh/read_mesh_files.F90 b/src/decompose_mesh/read_mesh_files.F90 index 0b7c3ab43..389be55e2 100644 --- a/src/decompose_mesh/read_mesh_files.F90 +++ b/src/decompose_mesh/read_mesh_files.F90 @@ -38,6 +38,8 @@ subroutine read_mesh_files() use fault_scotch, only: ANY_FAULT,read_fault_files,save_nodes_coords,close_faults + use constants, only: FNAME_WAVEFIELD_DISCONTINUITY_INTERFACE + implicit none ! local parameters @@ -46,7 +48,7 @@ subroutine read_mesh_files() ! poroelastic parameters read in a new file double precision :: rhos,rhof,phi,tort,kxx,kxy,kxz,kyy,kyz,kzz,kappas,kappaf,kappafr,eta,mufr integer(kind=8) :: nspec_long - integer :: inode,idummy + integer :: inode,idummy,ib integer :: ispec,ispec2D,ispec_CPML,ier character(len=MAX_STRING_LEN) :: line @@ -566,6 +568,32 @@ subroutine read_mesh_files() endif enddo + ! reads in wavefield discontinuity elements + if (IS_WAVEFIELD_DISCONTINUITY) then + open(unit=IIN_DB, file=localpath_name(1:len_trim(localpath_name))//& + '/'//trim(FNAME_WAVEFIELD_DISCONTINUITY_INTERFACE), & + status='old', form='formatted',iostat=ier) + if (ier /= 0) then + stop 'no wavefield discontinuity interface file' + endif + ier = 0 + nb_wd = 0 + do while (ier == 0) + read(IIN_DB, '(a)', iostat=ier) line + if (ier /= 0) exit + nb_wd = nb_wd + 1 + enddo + close(IIN_DB) + allocate(boundary_to_ispec_wd(nb_wd), side_wd(nb_wd)) + open(unit=IIN_DB, file=localpath_name(1:len_trim(localpath_name))//& + '/'//trim(FNAME_WAVEFIELD_DISCONTINUITY_INTERFACE), & + status='old', form='formatted',iostat=ier) + do ib = 1, nb_wd + read(IIN_DB, *, iostat=ier) boundary_to_ispec_wd(ib), side_wd(ib) + enddo + close(IIN_DB) + endif + ! reads in absorbing boundary files open(unit=IIN_DB, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_xmin', & status='old', form='formatted',iostat=ier) diff --git a/src/decompose_mesh/write_mesh_databases.F90 b/src/decompose_mesh/write_mesh_databases.F90 index f4745c690..7c42244a2 100644 --- a/src/decompose_mesh/write_mesh_databases.F90 +++ b/src/decompose_mesh/write_mesh_databases.F90 @@ -108,10 +108,12 @@ subroutine write_mesh_databases() ! writes out spectral element indices write(IIN_database) nspec_local + call write_partition_database(IIN_database, ipart, nspec_local, nspec, elmnts, & glob2loc_elmnts, glob2loc_nodes_nparts, & glob2loc_nodes_parts, glob2loc_nodes, part, mat, NGNOD, 2) + ! writes out absorbing/free-surface boundaries call write_boundaries_database(IIN_database, ipart, nspec, nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, & nspec2D_ymax, nspec2D_bottom, nspec2D_top, & @@ -199,6 +201,16 @@ subroutine write_mesh_databases() enddo endif + + !! setting up wavefield discontinuity boundary + if (IS_WAVEFIELD_DISCONTINUITY) then + do ipart = 0, nparts-1 + call write_wavefield_discontinuity_database(ipart, outputpath_name, & + nb_wd, boundary_to_ispec_wd, side_wd, & + nspec, glob2loc_elmnts, part) + enddo + endif + ! cleanup deallocate(CPML_to_spec,stat=ier); if (ier /= 0) stop 'Error deallocating array CPML_to_spec' deallocate(CPML_regions,stat=ier); if (ier /= 0) stop 'Error deallocating array CPML_regions' diff --git a/src/generate_databases/create_regions_mesh.f90 b/src/generate_databases/create_regions_mesh.f90 index e09ee313c..8286ea7c6 100644 --- a/src/generate_databases/create_regions_mesh.f90 +++ b/src/generate_databases/create_regions_mesh.f90 @@ -62,6 +62,12 @@ subroutine create_regions_mesh() nnodes_coords_open,nodes_coords_open,ANY_FAULT_IN_THIS_PROC, & ANY_FAULT + !! setup wavefield discontinuity interface + use shared_parameters, only: IS_WAVEFIELD_DISCONTINUITY + use wavefield_discontinuity_generate_databases, only: & + setup_boundary_wavefield_discontinuity, & + read_partition_files_wavefield_discontinuity + implicit none ! local parameters @@ -329,6 +335,18 @@ subroutine create_regions_mesh() ! user output call print_timing() + ! setup wavefield discontinuity interface + if (IS_WAVEFIELD_DISCONTINUITY) then + call synchronize_all() + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) ' ...setting up wavefield discontinuity boundary ' + call flush_IMAIN() + endif + call read_partition_files_wavefield_discontinuity() + call setup_boundary_wavefield_discontinuity() + endif + ! saves the binary mesh files call synchronize_all() if (myrank == 0) then diff --git a/src/generate_databases/rules.mk b/src/generate_databases/rules.mk index b413ddf91..1f4600aa6 100644 --- a/src/generate_databases/rules.mk +++ b/src/generate_databases/rules.mk @@ -45,6 +45,7 @@ generate_databases_OBJECTS = \ $O/generate_databases_par.gen_mod.o \ $O/calc_jacobian.gen.o \ $O/fault_generate_databases.gen.o \ + $O/wavefield_discontinuity_generate_databases.gen.o \ $O/create_mass_matrices.gen.o \ $O/create_regions_mesh.gen.o \ $O/finalize_databases.gen.o \ @@ -91,6 +92,7 @@ generate_databases_MODULES = \ $(FC_MODDIR)/model_sep_mod.$(FC_MODEXT) \ $(FC_MODDIR)/model_tomography_par.$(FC_MODEXT) \ $(FC_MODDIR)/salton_trough_par.$(FC_MODEXT) \ + $(FC_MODDIR)/wavefield_discontinuity_generate_databases.$(FC_MODEXT) \ $(EMPTY_MACRO) @@ -229,6 +231,10 @@ $O/lts_generate_databases.gen.o: $O/fault_generate_databases.gen.o ## kdtree & faults dependency $O/setup_mesh_adjacency.gen.o: $O/search_kdtree.shared.o $O/fault_generate_databases.gen.o +## wavefield discontinuity +$O/create_regions_mesh.gen.o: $O/wavefield_discontinuity_generate_databases.gen.o +$O/save_arrays_solver.gen.o: $O/wavefield_discontinuity_generate_databases.gen.o + ####################################### diff --git a/src/generate_databases/save_arrays_solver.F90 b/src/generate_databases/save_arrays_solver.F90 index 01b8cb10a..59cde021b 100644 --- a/src/generate_databases/save_arrays_solver.F90 +++ b/src/generate_databases/save_arrays_solver.F90 @@ -69,6 +69,11 @@ subroutine save_arrays_solver_mesh() use shared_parameters, only: ADIOS_FOR_MESH,HDF5_ENABLED + !! setup wavefield discontinuity interface + use shared_parameters, only: IS_WAVEFIELD_DISCONTINUITY + use wavefield_discontinuity_generate_databases, only: & + save_arrays_solver_mesh_wavefield_discontinuity + implicit none ! local parameters @@ -408,6 +413,11 @@ subroutine save_arrays_solver_mesh() call save_arrays_solver_injection_boundary() endif + !! setup wavefield discontinuity interface + if (IS_WAVEFIELD_DISCONTINUITY) then + call save_arrays_solver_mesh_wavefield_discontinuity() + endif + ! synchronizes processes call synchronize_all() diff --git a/src/generate_databases/wavefield_discontinuity_generate_databases.f90 b/src/generate_databases/wavefield_discontinuity_generate_databases.f90 new file mode 100644 index 000000000..12c3292d2 --- /dev/null +++ b/src/generate_databases/wavefield_discontinuity_generate_databases.f90 @@ -0,0 +1,693 @@ +module wavefield_discontinuity_generate_databases + use constants, only: CUSTOM_REAL + + !! boundary of wavefield discontinuity, read from database file + integer :: nb_wd + + !! boundary_to_ispec_wd(nb_wd) + !! the element the boundary belongs to, read from database file + !! each point on the boundary belongs to two sides of the boundary + !! here the element must be on the inner side of the boundary + integer, dimension(:), allocatable :: boundary_to_ispec_wd + + !! side_wd(nb_wd) + !! integers specifying which side the boundary is in the element + !! read from database file + !! side_wd = 1--8: only one vertex is on the boundary + !! side_wd = 9--20: only one edge is on the boundary + !! side_wd = 21--26: one face is on the boundary + integer, dimension(:), allocatable :: side_wd + + !! ispec_to_elem_wd(NSPEC_AB) + !! ispec_to_elem_wd(ispec) = ispec_wd (0 if element not belong to boundary) + !! written in solver database and used in solver + integer, dimension(:), allocatable :: ispec_to_elem_wd + + !! number of distinct gll points on the boundary + !! written in solver database and used in solver + integer :: nglob_wd + + !! number of elements on the inner side of the boundary + !! written in solver database and used in solver + integer :: nspec_wd + + !! ibool_wd(NGLLX, NGLLY, NGLLZ, nspec_wd) + !! ibool_wd(i,j,k,ispec_wd) = iglob_wd (0 if point not on boundary) + !! written in solver database and used in solver + integer, dimension(:,:,:,:), allocatable :: ibool_wd + + !! boundary_to_iglob_wd(nglob_wd) + !! boundary_to_iglob_wd(iglob_wd) = iglob + !! written in solver database and used in solver + integer, dimension(:), allocatable :: boundary_to_iglob_wd + + !! mass_in_wd(nglob_wd) + !! mass matrix on the inner side of the boundary + !! note that it is not assembled over processors + !! written in solver database and used in solver + real(kind=CUSTOM_REAL), dimension(:), allocatable :: mass_in_wd + + !! number of faces on the boundary + !! written in solver database and used in solver + integer :: nfaces_wd + + !! face_ijk_wd(NDIM, NGLLSQUARE, nfaces_wd) + !! written in solver database and used in solver + integer, dimension(:,:,:), allocatable :: face_ijk_wd + + !! face_ispec_wd(nfaces_wd) + !! written in solver database and used in solver + integer, dimension(:), allocatable :: face_ispec_wd + + !! face_normal_wd(NDIM, NGLLSQUARE, nfaces_wd) + !! written in solver database and used in solver + real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: face_normal_wd + + !! face_jacobian2Dw_wd(NGLLSQUARE, nfaces_wd) + !! written in solver database and used in solver + real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: face_jacobian2dw_wd + +contains + subroutine read_partition_files_wavefield_discontinuity() + use constants, only: IFILE_WAVEFIELD_DISCONTINUITY, & + FNAME_WAVEFIELD_DISCONTINUITY_MESH + use generate_databases_par, only: prname + implicit none + !integer :: ier + open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & + file=trim(prname)//trim(FNAME_WAVEFIELD_DISCONTINUITY_MESH), & + action='read', form='unformatted') + read(IFILE_WAVEFIELD_DISCONTINUITY) nb_wd + allocate(boundary_to_ispec_wd(nb_wd), side_wd(nb_wd)) + read(IFILE_WAVEFIELD_DISCONTINUITY) boundary_to_ispec_wd + read(IFILE_WAVEFIELD_DISCONTINUITY) side_wd + close(IFILE_WAVEFIELD_DISCONTINUITY) + end subroutine read_partition_files_wavefield_discontinuity + + subroutine save_arrays_solver_mesh_wavefield_discontinuity() + use constants, only: IFILE_WAVEFIELD_DISCONTINUITY, & + FNAME_WAVEFIELD_DISCONTINUITY_DATABASE + use generate_databases_par, only: prname + implicit none + open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & + file=trim(prname)//trim(FNAME_WAVEFIELD_DISCONTINUITY_DATABASE), & + action='write', form='unformatted') + write(IFILE_WAVEFIELD_DISCONTINUITY) ispec_to_elem_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) nglob_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) nspec_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) ibool_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) boundary_to_iglob_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) mass_in_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) nfaces_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) face_ijk_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) face_ispec_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) face_normal_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) face_jacobian2dw_wd + close(IFILE_WAVEFIELD_DISCONTINUITY) + deallocate(boundary_to_ispec_wd, side_wd) + deallocate(ispec_to_elem_wd, ibool_wd, boundary_to_iglob_wd, mass_in_wd, & + face_ijk_wd, face_ispec_wd, face_normal_wd, face_jacobian2dw_wd) + end subroutine save_arrays_solver_mesh_wavefield_discontinuity + + subroutine setup_boundary_wavefield_discontinuity() + use generate_databases_par, only: NDIM, NGLLX, NGLLY, NGLLZ, CUSTOM_REAL, & + NGLLSQUARE, ibool, NSPEC_AB + + implicit none + integer :: i1, i2, j1, j2, k1, k2, i, j, k + integer :: ib, iside, ispec, iglob, ispec_wd, iglob_wd + logical :: is_face + integer, dimension(:), allocatable :: elem_list_temp, gllp_list_temp + real(kind=CUSTOM_REAL) :: val_mass + allocate(ispec_to_elem_wd(NSPEC_AB)) + allocate(elem_list_temp(nb_wd), gllp_list_temp(nb_wd*NGLLSQUARE)) + nspec_wd = 0 + nglob_wd = 0 + nfaces_wd = 0 + ispec_to_elem_wd(:) = 0 + elem_list_temp(:) = 0 + gllp_list_temp(:) = 0 + do ib = 1, nb_wd + ispec = boundary_to_ispec_wd(ib) + iside = side_wd(ib) + call get_points_boundary_wd(iside, i1, i2, j1, j2, k1, k2, is_face) + call find_point_wd(ispec, elem_list_temp, nspec_wd, ispec_wd) + if (ispec_wd == 0) then + ispec_to_elem_wd(ispec) = nspec_wd + 1 + elem_list_temp(nspec_wd+1) = ispec + nspec_wd = nspec_wd + 1 + endif + do i=i1,i2; do j=j1,j2; do k=k1,k2 + iglob = ibool(i,j,k,ispec) + call find_point_wd(iglob, gllp_list_temp, nglob_wd, iglob_wd) + if (iglob_wd == 0) then + gllp_list_temp(nglob_wd+1) = iglob + nglob_wd = nglob_wd + 1 + endif + enddo; enddo; enddo + if (is_face) nfaces_wd = nfaces_wd + 1 + enddo + allocate(ibool_wd(NGLLX, NGLLY, NGLLZ, nspec_wd), & + boundary_to_iglob_wd(nglob_wd), & + mass_in_wd(nglob_wd), & + face_ijk_wd(NDIM, NGLLSQUARE, nfaces_wd), & + face_ispec_wd(nfaces_wd), & + face_normal_wd(NDIM, NGLLSQUARE, nfaces_wd), & + face_jacobian2Dw_wd(NGLLSQUARE, nfaces_wd)) + ibool_wd(:,:,:,:) = 0 + mass_in_wd(:) = 0.0 + nfaces_wd = 0 + do ib = 1, nb_wd + ispec = boundary_to_ispec_wd(ib) + iside = side_wd(ib) + call get_points_boundary_wd(iside, i1, i2, j1, j2, k1, k2, is_face) + call find_point_wd(ispec, elem_list_temp, nspec_wd, ispec_wd) + do i=i1,i2; do j=j1,j2; do k=k1,k2 + iglob = ibool(i,j,k,ispec) + call find_point_wd(iglob, gllp_list_temp, nglob_wd, iglob_wd) + ibool_wd(i,j,k,ispec_wd) = iglob_wd + boundary_to_iglob_wd(iglob_wd) = iglob + !call get_mass_wd(i, j, k, ispec, val_mass) + !mass_in_wd(iglob_wd) = mass_in_wd(iglob_wd) + val_mass + enddo; enddo; enddo + if (is_face) then + nfaces_wd = nfaces_wd + 1 + face_ispec_wd(nfaces_wd) = ispec + call get_face_wd(ispec, iside, face_ijk_wd(:,:,nfaces_wd), & + face_normal_wd(:,:,nfaces_wd), & + face_jacobian2Dw_wd(:,nfaces_wd)) + endif + enddo + do ispec = 1, NSPEC_AB + ispec_wd = ispec_to_elem_wd(ispec) + if (ispec_wd > 0) then + do k=1,NGLLZ; do j=1,NGLLY; do i=1, NGLLX + iglob_wd= ibool_wd(i,j,k,ispec_wd) + if (iglob_wd > 0) then + call get_mass_wd(i, j, k, ispec, val_mass) + mass_in_wd(iglob_wd) = mass_in_wd(iglob_wd) + val_mass + endif + enddo; enddo; enddo + endif + enddo + call write_discontinuity_surface_file() + end subroutine setup_boundary_wavefield_discontinuity + + subroutine write_discontinuity_surface_file() + use constants, only: NGLLX,NGLLY,NGLLZ,NDIM,NGLLSQUARE,CUSTOM_REAL,& + IFILE_WAVEFIELD_DISCONTINUITY + use generate_databases_par, only: prname, myrank, LOCAL_PATH, & + nodes_coords_ext_mesh, elmnts_ext_mesh, NGNOD + use create_regions_mesh_ext_par, only: xstore_unique, ystore_unique, & + zstore_unique, & + xigll,yigll,zigll + implicit none + double precision, parameter :: ONE_SHRINK = 0.999 + integer :: i,j,k,ispec,iglob,igll,iglob_wd,iface_wd, ia + real(CUSTOM_REAL) :: tx, ty, tz + double precision :: shape3D_shrink(NGNOD,NGLLX,NGLLY,NGLLZ) + double precision :: dershape3D_shrink(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ) + double precision :: xigll_shrink(NGLLX), yigll_shrink(NGLLY), & + zigll_shrink(NGLLZ) + double precision :: xstore_shrink(NGLLX,NGLLY,NGLLZ), & + ystore_shrink(NGLLX,NGLLY,NGLLZ), & + zstore_shrink(NGLLX,NGLLY,NGLLZ) + double precision, dimension(NGNOD) :: xelm,yelm,zelm + call create_name_database(prname,myrank,LOCAL_PATH) + open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & + file=prname(1:len_trim(prname))//'wavefield_discontinuity_points',& + action='write', form='formatted') + do iglob_wd = 1, nglob_wd + iglob = boundary_to_iglob_wd(iglob_wd) + write(IFILE_WAVEFIELD_DISCONTINUITY, '(4e20.5)') xstore_unique(iglob), & + ystore_unique(iglob), zstore_unique(iglob), mass_in_wd(iglob_wd) + enddo + close(IFILE_WAVEFIELD_DISCONTINUITY) + xigll_shrink(:) = xigll(:) + xigll_shrink(1) = xigll_shrink(1) * ONE_SHRINK + xigll_shrink(NGLLX) = xigll_shrink(NGLLX) * ONE_SHRINK + yigll_shrink(:) = yigll(:) + yigll_shrink(1) = yigll_shrink(1) * ONE_SHRINK + yigll_shrink(NGLLX) = yigll_shrink(NGLLX) * ONE_SHRINK + zigll_shrink(:) = zigll(:) + zigll_shrink(1) = zigll_shrink(1) * ONE_SHRINK + zigll_shrink(NGLLX) = zigll_shrink(NGLLX) * ONE_SHRINK + call get_shape3D(shape3D_shrink, dershape3D_shrink, & + xigll_shrink, yigll_shrink, zigll_shrink, NGNOD, NGLLX, NGLLY, NGLLZ) + open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & + file=prname(1:len_trim(prname))//'wavefield_discontinuity_faces',& + action='write', form='formatted') + do iface_wd = 1, nfaces_wd + ispec = face_ispec_wd(iface_wd) + do ia = 1,NGNOD + iglob = elmnts_ext_mesh(ia,ispec) + xelm(ia) = nodes_coords_ext_mesh(1,iglob) + yelm(ia) = nodes_coords_ext_mesh(2,iglob) + zelm(ia) = nodes_coords_ext_mesh(3,iglob) + enddo + call calc_coords(xstore_shrink(1,1,1), ystore_shrink(1,1,1),& + zstore_shrink(1,1,1), xelm,yelm,zelm,shape3D_shrink) + do igll = 1, NGLLSQUARE + i = face_ijk_wd(1,igll,iface_wd) + j = face_ijk_wd(2,igll,iface_wd) + k = face_ijk_wd(3,igll,iface_wd) + !iglob = ibool(i,j,k,ispec) + tx = face_normal_wd(1,igll,iface_wd) + ty = face_normal_wd(2,igll,iface_wd) + tz = face_normal_wd(3,igll,iface_wd) + !write(IFILE_WAVEFIELD_DISCONTINUITY, '(7e20.5)') xstore_unique(iglob), & + ! ystore_unique(iglob), zstore_unique(iglob), tx, ty, tz, & + ! face_jacobian2Dw_wd(igll, iface_wd) + write(IFILE_WAVEFIELD_DISCONTINUITY, '(7e20.5)') xstore_shrink(i,j,k), & + ystore_shrink(i,j,k), zstore_shrink(i,j,k), tx, ty, tz, & + face_jacobian2Dw_wd(igll, iface_wd) + enddo + enddo + end subroutine write_discontinuity_surface_file + + subroutine get_face_wd(ispec, iside, face_ijk_b, face_normal_b, & + face_jacobian2Dw_b) + use constants, only: NGLLX,NGLLY,NGLLZ,NDIM,NGNOD2D_FOUR_CORNERS, & + CUSTOM_REAL, NGLLSQUARE + use generate_databases_par, only: NGNOD2D, ibool, NSPEC_AB, NGLOB_AB, & + xstore, ystore, zstore + use create_regions_mesh_ext_par, only: xstore_unique, ystore_unique, & + zstore_unique, nglob_unique, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz + implicit none + integer, intent(in) :: ispec, iside + integer :: iface, igll, i, j + integer :: ijk_face(NDIM, NGLLX, NGLLY) + real(kind=CUSTOM_REAL) :: normal_face(NDIM, NGLLX, NGLLY), & + jacobian2Dw_face(NGLLX, NGLLY) + real(kind=CUSTOM_REAL), dimension(NDIM) :: lnormal + integer, intent(out) :: face_ijk_b(NDIM, NGLLSQUARE) + real(kind=CUSTOM_REAL), intent(out) :: face_normal_b(NDIM, NGLLSQUARE), & + face_jacobian2Dw_b(NGLLSQUARE) + real(kind=CUSTOM_REAL), dimension(NGNOD2D_FOUR_CORNERS) :: & + xcoord,ycoord,zcoord + if (nglob_unique /= NGLOB_AB) & + call exit_MPI_without_rank('nglob_unique not equal to NGLOB_AB') + select case (iside) + case (21) + ! bottom + xcoord(1) = xstore(1,1,1,ispec) + xcoord(2) = xstore(1,NGLLY,1,ispec) + xcoord(3) = xstore(NGLLX,NGLLY,1,ispec) + xcoord(4) = xstore(NGLLX,1,1,ispec) + ycoord(1) = ystore(1,1,1,ispec) + ycoord(2) = ystore(1,NGLLY,1,ispec) + ycoord(3) = ystore(NGLLX,NGLLY,1,ispec) + ycoord(4) = ystore(NGLLX,1,1,ispec) + zcoord(1) = zstore(1,1,1,ispec) + zcoord(2) = zstore(1,NGLLY,1,ispec) + zcoord(3) = zstore(NGLLX,NGLLY,1,ispec) + zcoord(4) = zstore(NGLLX,1,1,ispec) + ! sets face id of reference element associated with this face + call get_element_face_id(ispec,xcoord,ycoord,zcoord, & + ibool,NSPEC_AB,NGLOB_AB, & + xstore_unique,ystore_unique,zstore_unique,iface) + if (iface /= 5) & + call exit_MPI_without_rank('incorrect iface') + ! ijk indices of GLL points on face + call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLY) + ! weighted jacobian and normal + call get_jacobian_boundary_face(NSPEC_AB, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) + ! normal convention: points away from element + ! switch normal direction if necessary + do j = 1,NGLLY + do i = 1,NGLLX + lnormal(:) = normal_face(:,i,j) + call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & + ibool,NSPEC_AB,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & + lnormal) + normal_face(:,i,j) = lnormal(:) + enddo + enddo + ! GLL points -- assuming NGLLX = NGLLY = NGLLZ + igll = 0 + do j = 1,NGLLY + do i = 1,NGLLX + igll = igll+1 + face_ijk_b(:,igll) = ijk_face(:,i,j) + face_jacobian2Dw_b(igll) = jacobian2Dw_face(i,j) + face_normal_b(:,igll) = normal_face(:,i,j) + enddo + enddo + case (25) + ! left, xmin + xcoord(1) = xstore(1,1,1,ispec) + xcoord(2) = xstore(1,1,NGLLZ,ispec) + xcoord(3) = xstore(1,NGLLY,NGLLZ,ispec) + xcoord(4) = xstore(1,NGLLY,1,ispec) + ycoord(1) = ystore(1,1,1,ispec) + ycoord(2) = ystore(1,1,NGLLZ,ispec) + ycoord(3) = ystore(1,NGLLY,NGLLZ,ispec) + ycoord(4) = ystore(1,NGLLY,1,ispec) + zcoord(1) = zstore(1,1,1,ispec) + zcoord(2) = zstore(1,1,NGLLZ,ispec) + zcoord(3) = zstore(1,NGLLY,NGLLZ,ispec) + zcoord(4) = zstore(1,NGLLY,1,ispec) + ! sets face id of reference element associated with this face + call get_element_face_id(ispec,xcoord,ycoord,zcoord, & + ibool,NSPEC_AB,NGLOB_AB, & + xstore_unique,ystore_unique,zstore_unique,iface) + if (iface /= 1) & + call exit_MPI_without_rank('incorrect iface') + ! ijk indices of GLL points on face + call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLZ) + ! weighted jacobian and normal + call get_jacobian_boundary_face(NSPEC_AB, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) + ! normal convention: points away from element + ! switch normal direction if necessary + do j = 1,NGLLZ + do i = 1,NGLLX + lnormal(:) = normal_face(:,i,j) + call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & + ibool,NSPEC_AB,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & + lnormal) + normal_face(:,i,j) = lnormal(:) + enddo + enddo + ! GLL points -- assuming NGLLX = NGLLY = NGLLZ + igll = 0 + do j = 1,NGLLZ + do i = 1,NGLLX + igll = igll+1 + face_ijk_b(:,igll) = ijk_face(:,i,j) + face_jacobian2Dw_b(igll) = jacobian2Dw_face(i,j) + face_normal_b(:,igll) = normal_face(:,i,j) + enddo + enddo + case (22) + ! front, ymin + xcoord(1) = xstore(1,1,1,ispec) + xcoord(2) = xstore(NGLLX,1,1,ispec) + xcoord(3) = xstore(NGLLX,1,NGLLZ,ispec) + xcoord(4) = xstore(1,1,NGLLZ,ispec) + ycoord(1) = ystore(1,1,1,ispec) + ycoord(2) = ystore(NGLLX,1,1,ispec) + ycoord(3) = ystore(NGLLX,1,NGLLZ,ispec) + ycoord(4) = ystore(1,1,NGLLZ,ispec) + zcoord(1) = zstore(1,1,1,ispec) + zcoord(2) = zstore(NGLLX,1,1,ispec) + zcoord(3) = zstore(NGLLX,1,NGLLZ,ispec) + zcoord(4) = zstore(1,1,NGLLZ,ispec) + ! sets face id of reference element associated with this face + call get_element_face_id(ispec,xcoord,ycoord,zcoord, & + ibool,NSPEC_AB,NGLOB_AB, & + xstore_unique,ystore_unique,zstore_unique,iface) + if (iface /= 3) & + call exit_MPI_without_rank('incorrect iface') + ! ijk indices of GLL points on face + call get_element_face_gll_indices(iface,ijk_face,NGLLY,NGLLZ) + ! weighted jacobian and normal + call get_jacobian_boundary_face(NSPEC_AB, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) + ! normal convention: points away from element + ! switch normal direction if necessary + do j = 1,NGLLZ + do i = 1,NGLLY + lnormal(:) = normal_face(:,i,j) + call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & + ibool,NSPEC_AB,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & + lnormal) + normal_face(:,i,j) = lnormal(:) + enddo + enddo + ! GLL points -- assuming NGLLX = NGLLY = NGLLZ + igll = 0 + do j = 1,NGLLZ + do i = 1,NGLLY + igll = igll+1 + face_ijk_b(:,igll) = ijk_face(:,i,j) + face_jacobian2Dw_b(igll) = jacobian2Dw_face(i,j) + face_normal_b(:,igll) = normal_face(:,i,j) + enddo + enddo + case (23) + ! right, xmax + xcoord(1) = xstore(NGLLX,1,1,ispec) + xcoord(2) = xstore(NGLLX,NGLLY,1,ispec) + xcoord(3) = xstore(NGLLX,NGLLY,NGLLZ,ispec) + xcoord(4) = xstore(NGLLX,1,NGLLZ,ispec) + ycoord(1) = ystore(NGLLX,1,1,ispec) + ycoord(2) = ystore(NGLLX,NGLLY,1,ispec) + ycoord(3) = ystore(NGLLX,NGLLY,NGLLZ,ispec) + ycoord(4) = ystore(NGLLX,1,NGLLZ,ispec) + zcoord(1) = zstore(NGLLX,1,1,ispec) + zcoord(2) = zstore(NGLLX,NGLLY,1,ispec) + zcoord(3) = zstore(NGLLX,NGLLY,NGLLZ,ispec) + zcoord(4) = zstore(NGLLX,1,NGLLZ,ispec) + ! sets face id of reference element associated with this face + call get_element_face_id(ispec,xcoord,ycoord,zcoord, & + ibool,NSPEC_AB,NGLOB_AB, & + xstore_unique,ystore_unique,zstore_unique,iface) + if (iface /= 2) & + call exit_MPI_without_rank('incorrect iface') + ! ijk indices of GLL points on face + call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLZ) + ! weighted jacobian and normal + call get_jacobian_boundary_face(NSPEC_AB, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) + ! normal convention: points away from element + ! switch normal direction if necessary + do j = 1,NGLLZ + do i = 1,NGLLX + lnormal(:) = normal_face(:,i,j) + call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & + ibool,NSPEC_AB,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & + lnormal) + normal_face(:,i,j) = lnormal(:) + enddo + enddo + ! GLL points -- assuming NGLLX = NGLLY = NGLLZ + igll = 0 + do j = 1,NGLLZ + do i = 1,NGLLX + igll = igll+1 + face_ijk_b(:,igll) = ijk_face(:,i,j) + face_jacobian2Dw_b(igll) = jacobian2Dw_face(i,j) + face_normal_b(:,igll) = normal_face(:,i,j) + enddo + enddo + case (24) + ! back, ymax + xcoord(1) = xstore(1,NGLLY,1,ispec) + xcoord(2) = xstore(1,NGLLY,NGLLZ,ispec) + xcoord(3) = xstore(NGLLX,NGLLY,NGLLZ,ispec) + xcoord(4) = xstore(NGLLX,NGLLY,1,ispec) + ycoord(1) = ystore(1,NGLLY,1,ispec) + ycoord(2) = ystore(1,NGLLY,NGLLZ,ispec) + ycoord(3) = ystore(NGLLX,NGLLY,NGLLZ,ispec) + ycoord(4) = ystore(NGLLX,NGLLY,1,ispec) + zcoord(1) = zstore(1,NGLLY,1,ispec) + zcoord(2) = zstore(1,NGLLY,NGLLZ,ispec) + zcoord(3) = zstore(NGLLX,NGLLY,NGLLZ,ispec) + zcoord(4) = zstore(NGLLX,NGLLY,1,ispec) + ! sets face id of reference element associated with this face + call get_element_face_id(ispec,xcoord,ycoord,zcoord, & + ibool,NSPEC_AB,NGLOB_AB, & + xstore_unique,ystore_unique,zstore_unique,iface) + if (iface /= 4) & + call exit_MPI_without_rank('incorrect iface') + ! ijk indices of GLL points on face + call get_element_face_gll_indices(iface,ijk_face,NGLLY,NGLLZ) + ! weighted jacobian and normal + call get_jacobian_boundary_face(NSPEC_AB, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) + ! normal convention: points away from element + ! switch normal direction if necessary + do j = 1,NGLLZ + do i = 1,NGLLY + lnormal(:) = normal_face(:,i,j) + call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & + ibool,NSPEC_AB,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & + lnormal) + normal_face(:,i,j) = lnormal(:) + enddo + enddo + ! GLL points -- assuming NGLLX = NGLLY = NGLLZ + igll = 0 + do j = 1,NGLLZ + do i = 1,NGLLY + igll = igll+1 + face_ijk_b(:,igll) = ijk_face(:,i,j) + face_jacobian2Dw_b(igll) = jacobian2Dw_face(i,j) + face_normal_b(:,igll) = normal_face(:,i,j) + enddo + enddo + case (26) + ! top, zmax + xcoord(1) = xstore(1,1,NGLLZ,ispec) + xcoord(2) = xstore(NGLLX,1,NGLLZ,ispec) + xcoord(3) = xstore(NGLLX,NGLLY,NGLLZ,ispec) + xcoord(4) = xstore(1,NGLLY,NGLLZ,ispec) + ycoord(1) = ystore(1,1,NGLLZ,ispec) + ycoord(2) = ystore(NGLLX,1,NGLLZ,ispec) + ycoord(3) = ystore(NGLLX,NGLLY,NGLLZ,ispec) + ycoord(4) = ystore(1,NGLLY,NGLLZ,ispec) + zcoord(1) = zstore(1,1,NGLLZ,ispec) + zcoord(2) = zstore(NGLLX,1,NGLLZ,ispec) + zcoord(3) = zstore(NGLLX,NGLLY,NGLLZ,ispec) + zcoord(4) = zstore(1,NGLLY,NGLLZ,ispec) + ! sets face id of reference element associated with this face + call get_element_face_id(ispec,xcoord,ycoord,zcoord, & + ibool,NSPEC_AB,NGLOB_AB, & + xstore_unique,ystore_unique,zstore_unique,iface) + if (iface /= 6) & + call exit_MPI_without_rank('incorrect iface') + ! ijk indices of GLL points on face + call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLY) + ! weighted jacobian and normal + call get_jacobian_boundary_face(NSPEC_AB, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) + ! normal convention: points away from element + ! switch normal direction if necessary + do j = 1,NGLLY + do i = 1,NGLLX + lnormal(:) = normal_face(:,i,j) + call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & + ibool,NSPEC_AB,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & + lnormal) + normal_face(:,i,j) = lnormal(:) + enddo + enddo + ! GLL points -- assuming NGLLX = NGLLY = NGLLZ + igll = 0 + do j = 1,NGLLY + do i = 1,NGLLX + igll = igll+1 + face_ijk_b(:,igll) = ijk_face(:,i,j) + face_jacobian2Dw_b(igll) = jacobian2Dw_face(i,j) + face_normal_b(:,igll) = normal_face(:,i,j) + enddo + enddo + end select + end subroutine get_face_wd + + subroutine get_mass_wd(i, j, k, ispec, val_mass) + use generate_databases_par, only: CUSTOM_REAL + use create_regions_mesh_ext_par, only: irregular_element_number, & + wxgll, wygll, wzgll, jacobianstore, & + rhostore, jacobian_regular + implicit none + integer :: i, j, k, ispec, ispec_irreg + real(kind=CUSTOM_REAL) :: val_mass, jacobianl + double precision :: weight + ispec_irreg = irregular_element_number(ispec) + if (ispec_irreg == 0) then + jacobianl = jacobian_regular + else + jacobianl = jacobianstore(i,j,k,ispec_irreg) + endif + weight = wxgll(i)*wygll(j)*wzgll(k) + val_mass = real(dble(jacobianl) * weight * dble(rhostore(i,j,k,ispec)), & + kind=CUSTOM_REAL) + end subroutine get_mass_wd + + subroutine find_point_wd(p, plist, len, ip) + implicit none + integer, intent(in) :: p, len + integer, dimension(len), intent(in) :: plist + integer, intent(out) :: ip + integer :: i + ip = 0 + do i = 1, len + if (plist(i) == p) then + ip = i + exit + endif + enddo + end subroutine find_point_wd + + subroutine get_points_boundary_wd(iside, i1, i2, j1, j2, k1, k2, is_face) + use generate_databases_par, only: NGLLX, NGLLY, NGLLZ + implicit none + integer, intent(in) :: iside + integer, intent(out) :: i1, i2, j1, j2, k1, k2 + logical, intent(out) :: is_face + integer, parameter :: NGNOD = 27, THREE = 3 + integer, dimension(NGNOD) :: iaddx, iaddy, iaddz + integer :: i, j, k + call usual_hex_nodes(NGNOD, iaddx, iaddy, iaddz) + iaddx = iaddx * 2 + 1 + iaddy = iaddy * 2 + 1 + iaddz = iaddz * 2 + 1 + i = iaddx(iside) + j = iaddy(iside) + k = iaddz(iside) + is_face = .false. + if ((iside >= 1) .and. (iside <= 8)) then + ! a vertex is on the boundary + i1 = i; i2 = i + j1 = j; j2 = j + k1 = k; k2 = k + else if ((iside >= 9) .and. (iside <= 20)) then + ! an edge is on the boundary + if (i == THREE) then + i1 = 1; i2 = NGLLX + j1 = j; j2 = j + k1 = k; k2 = k + else if (j == THREE) then + i1 = i; i2 = i + j1 = 1; j2 = NGLLY + k1 = k; k2 = k + else if (k == THREE) then + i1 = i; i2 = i + j1 = j; j2 = j + k1 = 1; k2 = NGLLZ + else + call exit_MPI_without_rank('incorrect wavefield discontinuity point') + endif + else if ((iside >= 21) .and. (iside <= 26)) then + ! a face is on the boundary + is_face = .true. + if (i /= THREE) then + i1 = i; i2 = i + j1 = 1; j2 = NGLLY + k1 = 1; k2 = NGLLZ + else if (j /= THREE) then + i1 = 1; i2 = NGLLX + j1 = j; j2 = j + k1 = 1; k2 = NGLLZ + else if (k /= THREE) then + i1 = 1; i2 = NGLLX + j1 = 1; j2 = NGLLY + k1 = k; k2 = k + else + call exit_MPI_without_rank('incorrect wavefield discontinuity point') + endif + else + call exit_MPI_without_rank('incorrect wavefield discontinuity point') + endif + end subroutine get_points_boundary_wd + +end module wavefield_discontinuity_generate_databases diff --git a/src/inverse_problem_for_model/rules.mk b/src/inverse_problem_for_model/rules.mk index c6de580b1..55e2322aa 100644 --- a/src/inverse_problem_for_model/rules.mk +++ b/src/inverse_problem_for_model/rules.mk @@ -111,6 +111,7 @@ inverse_problem_for_model_OBJECTS = \ inverse_problem_for_model_OBJECTS += \ $O/specfem3D_par.spec_module.o \ $O/asdf_data.spec_module.o \ + $O/wavefield_discontinuity_solver.spec.o \ $O/assemble_MPI_vector.spec.o \ $O/calendar.spec.o \ $O/check_stability.spec.o \ diff --git a/src/meshfem3D/create_meshfem_mesh.f90 b/src/meshfem3D/create_meshfem_mesh.f90 index 5162c8127..5fe925c07 100644 --- a/src/meshfem3D/create_meshfem_mesh.f90 +++ b/src/meshfem3D/create_meshfem_mesh.f90 @@ -82,6 +82,9 @@ subroutine create_meshfem_mesh() ! HDF5 file i/o use shared_parameters, only: HDF5_ENABLED + + !! setting up wavefield discontinuity interface + use shared_parameters, only: IS_WAVEFIELD_DISCONTINUITY implicit none @@ -143,6 +146,11 @@ subroutine create_meshfem_mesh() nspec,nglob, & prname,nodes_coords,ibool,ispec_material_id) + !! setting up wavefield discontinuity interface + if (IS_WAVEFIELD_DISCONTINUITY) then + call find_wavefield_discontinuity_elements() + endif + ! stores boundary informations call store_boundaries(iboun,nspec, & ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, & diff --git a/src/meshfem3D/get_wavefield_discontinuity.f90 b/src/meshfem3D/get_wavefield_discontinuity.f90 new file mode 100644 index 000000000..7f3f26090 --- /dev/null +++ b/src/meshfem3D/get_wavefield_discontinuity.f90 @@ -0,0 +1,512 @@ +!===================================================================== +! +! S p e c f e m 3 D +! ----------------- +! +! Main historical authors: Dimitri Komatitsch and Jeroen Tromp +! CNRS, France +! and Princeton University, USA +! (there are currently many more authors!) +! (c) October 2017 +! +! 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 3 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 write_wavefield_discontinuity_database() +! write data related to wavefield discontinuity to proc*_Database + use meshfem_par, only: nb_wd, boundary_to_ispec_wd, side_wd, prname + use constants, only: FNAME_WAVEFIELD_DISCONTINUITY_MESH, & + IFILE_WAVEFIELD_DISCONTINUITY + implicit none + ! local variables + open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & + file=prname(1:len_trim(prname))//& + trim(FNAME_WAVEFIELD_DISCONTINUITY_MESH), & + form='unformatted', action='write') + write(IFILE_WAVEFIELD_DISCONTINUITY) nb_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) boundary_to_ispec_wd + write(IFILE_WAVEFIELD_DISCONTINUITY) side_wd + close(IFILE_WAVEFIELD_DISCONTINUITY) + end subroutine write_wavefield_discontinuity_database + +! +!----------------------------------------------------------- +! + + subroutine write_wavefield_discontinuity_file() +! write wavefield discontinuity interfaces to an ascii file +! works only when NPROC = 1 + use meshfem_par, only: nb_wd, boundary_to_ispec_wd, side_wd + use constants, only: IFILE_WAVEFIELD_DISCONTINUITY, & + FNAME_WAVEFIELD_DISCONTINUITY_INTERFACE + implicit none + ! local variables + integer :: i + open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & + file='MESH/'//trim(FNAME_WAVEFIELD_DISCONTINUITY_INTERFACE), & + form='formatted', action='write') + do i = 1, nb_wd + write(IFILE_WAVEFIELD_DISCONTINUITY, *) boundary_to_ispec_wd(i), side_wd(i) + enddo + close(IFILE_WAVEFIELD_DISCONTINUITY) + end subroutine write_wavefield_discontinuity_file + +! +!----------------------------------------------------------- +! + + subroutine find_wavefield_discontinuity_elements() +! read the wavefield_discontinuity_box file +! find the wavefield discontinuity interface + use constants, only: CUSTOM_REAL + use constants, only: IFILE_WAVEFIELD_DISCONTINUITY, & + FNAME_WAVEFIELD_DISCONTINUITY_BOX + use meshfem_par, only: xstore,ystore,zstore,nspec + use meshfem_par, only: nb_wd, boundary_to_ispec_wd, side_wd + implicit none + ! local variables + integer :: boundary_to_ispec_wd_temp(6*nspec), side_wd_temp(6*nspec) + integer :: ispec, iside, i, j, k + logical :: is_boundary_wd, IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE + logical :: covered(26) + double precision :: x_min, x_max, y_min, y_max, z_min, z_max + double precision :: dx, dy, dz, x_mid, y_mid, z_mid, ratio_small=1.0e-6 + open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & + file=trim(FNAME_WAVEFIELD_DISCONTINUITY_BOX), & + form='formatted', action='read') + read(IFILE_WAVEFIELD_DISCONTINUITY, *) IS_TOP_WAVEFIELD_DISCONTINUITY + read(IFILE_WAVEFIELD_DISCONTINUITY, *) IS_EXTRAPOLATION_MODE + read(IFILE_WAVEFIELD_DISCONTINUITY, *) x_min + read(IFILE_WAVEFIELD_DISCONTINUITY, *) x_max + read(IFILE_WAVEFIELD_DISCONTINUITY, *) y_min + read(IFILE_WAVEFIELD_DISCONTINUITY, *) y_max + read(IFILE_WAVEFIELD_DISCONTINUITY, *) z_min + read(IFILE_WAVEFIELD_DISCONTINUITY, *) z_max + close(IFILE_WAVEFIELD_DISCONTINUITY) + nb_wd = 0 + do ispec = 1, nspec + covered(:) = .false. + x_mid = xstore(2,2,2,ispec) + y_mid = ystore(2,2,2,ispec) + z_mid = zstore(2,2,2,ispec) + dx = ratio_small * abs(xstore(3,1,1,ispec) - xstore(1,1,1,ispec)) + dy = ratio_small * abs(ystore(1,3,1,ispec) - ystore(1,1,1,ispec)) + dz = ratio_small * abs(zstore(1,1,3,ispec) - zstore(1,1,1,ispec)) + !! bottom + iside = 21; i = 2; j = 2; k = 1 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/9,10,11,12,1,2,3,4,21/)) = .true. + endif + endif + !! front + iside = 22; i = 2; j = 1; k = 2 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/9,13,14,17,1,2,5,6,22/)) = .true. + endif + endif + !! right + iside = 23; i = 3; j = 2; k = 2 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/10,14,15,18,2,3,6,7,23/)) = .true. + endif + endif + !! back + iside = 24; i = 2; j = 3; k = 2 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/11,15,16,19,3,4,7,8,24/)) = .true. + endif + endif + !! left + iside = 25; i = 1; j = 2; k = 2 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/12,13,16,20,1,4,5,8,25/)) = .true. + endif + endif + !! top + iside = 26; i = 2; j = 2; k = 3 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/17,18,19,20,5,6,7,8,26/)) = .true. + endif + endif + !! front - bottom + iside = 9; i = 2; j = 1; k = 1 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/1,2,9/)) = .true. + endif + endif + !! right - bottom + iside = 10; i = 3; j = 2; k = 1 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/2,3,10/)) = .true. + endif + endif + !! back - bottom + iside = 11; i = 2; j = 3; k = 1 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/3,4,11/)) = .true. + endif + endif + !! left - bottom + iside = 12; i = 1; j = 2; k = 1 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/1,4,12/)) = .true. + endif + endif + !! left - front + iside = 13; i = 1; j = 1; k = 2 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/1,5,13/)) = .true. + endif + endif + !! right - front + iside = 14; i = 3; j = 1; k = 2 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/2,6,14/)) = .true. + endif + endif + !! right - back + iside = 15; i = 3; j = 3; k = 2 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/3,7,15/)) = .true. + endif + endif + !! left - front + iside = 16; i = 1; j = 3; k = 2 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/4,8,16/)) = .true. + endif + endif + !! front - top + iside = 17; i = 2; j = 1; k = 3 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/5,6,17/)) = .true. + endif + endif + !! right - top + iside = 18; i = 3; j = 2; k = 3 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/6,7,18/)) = .true. + endif + endif + !! back - top + iside = 19; i = 2; j = 3; k = 3 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/7,8,19/)) = .true. + endif + endif + !! left - bottom + iside = 20; i = 1; j = 2; k = 3 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered((/5,8,20/)) = .true. + endif + endif + !! left - front - bottom + iside = 1; i = 1; j = 1; k = 1 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered(1) = .true. + endif + endif + !! right - front - bottom + iside = 2; i = 3; j = 1; k = 1 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered(2) = .true. + endif + endif + !! right - back - bottom + iside = 3; i = 3; j = 3; k = 1 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered(3) = .true. + endif + endif + !! left - front - bottom + iside = 4; i = 1; j = 3; k = 1 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered(4) = .true. + endif + endif + !! left - front - top + iside = 5; i = 1; j = 1; k = 3 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered(5) = .true. + endif + endif + !! right - front - top + iside = 6; i = 3; j = 1; k = 3 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered(6) = .true. + endif + endif + !! right - back - top + iside = 7; i = 3; j = 3; k = 3 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered(7) = .true. + endif + endif + !! left - front - top + iside = 8; i = 1; j = 3; k = 3 + if (.not. covered(iside)) then + if (is_boundary_wd(xstore(i,j,k,ispec), ystore(i,j,k,ispec), & + zstore(i,j,k,ispec), x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE)) then + nb_wd = nb_wd + 1 + boundary_to_ispec_wd_temp(nb_wd) = ispec + side_wd_temp(nb_wd) = iside + covered(8) = .true. + endif + endif + enddo + allocate(boundary_to_ispec_wd(nb_wd), side_wd(nb_wd)) + boundary_to_ispec_wd(1:nb_wd) = boundary_to_ispec_wd_temp(1:nb_wd) + side_wd(1:nb_wd) = side_wd_temp(1:nb_wd) + end subroutine find_wavefield_discontinuity_elements + +! +!----------------------------------------------------------- +! + + logical function is_boundary_wd(x, y, z, x_min, x_max, y_min, y_max, & + z_min, z_max, x_mid, y_mid, z_mid, dx, dy, dz, & + IS_TOP_WAVEFIELD_DISCONTINUITY, & + IS_EXTRAPOLATION_MODE) + implicit none + double precision :: x_min, x_max, y_min, y_max, z_min, z_max + double precision :: x, y, z, x_mid, y_mid, z_mid, dx, dy, dz + logical :: IS_TOP_WAVEFIELD_DISCONTINUITY, IS_EXTRAPOLATION_MODE + is_boundary_wd = .false. + if (IS_EXTRAPOLATION_MODE) then + if (((x > x_min - dx) .and. (x < x_max + dx) .and. & + (y > y_min - dy) .and. (y < y_max + dy) .and. & + (z > z_min - dz) .and. (z < z_max + dz)) .and. & + ((x_mid < x_min) .or. (x_mid > x_max) .or. & + (y_mid < y_min) .or. (y_mid > y_max) .or. & + (z_mid < z_min) .or. (z_mid > z_max))) & + is_boundary_wd = .true. + else + if (((x < x_min + dx) .or. (x > x_max - dx) .or. & + (y < y_min + dy) .or. (y > y_max - dy) .or. & + (z < z_min + dz) .or. ((z > z_max - dz) .and. & + IS_TOP_WAVEFIELD_DISCONTINUITY)) .and. & + (x_mid > x_min) .and. (x_mid < x_max) .and. & + (y_mid > y_min) .and. (y_mid < y_max) .and. & + (z_mid > z_min) .and. ((.not. IS_TOP_WAVEFIELD_DISCONTINUITY) .or. & + (z_mid < z_max))) & + is_boundary_wd = .true. + endif + end function is_boundary_wd diff --git a/src/meshfem3D/meshfem3D_par.f90 b/src/meshfem3D/meshfem3D_par.f90 index 7c9049136..8faad973d 100644 --- a/src/meshfem3D/meshfem3D_par.f90 +++ b/src/meshfem3D/meshfem3D_par.f90 @@ -141,5 +141,22 @@ module meshfem_par ! name of the database file character(len=MAX_STRING_LEN) :: prname + !! boundary of wavefield discontinuity, read from database file + integer :: nb_wd + + !! boundary_to_ispec_wd(nb_wd) + !! the element the boundary belongs to, read from database file + !! each point on the boundary belongs to two sides of the boundary + !! here the element must be on the inner side of the boundary + integer, dimension(:), allocatable :: boundary_to_ispec_wd + + !! side_wd(nb_wd) + !! integers specifying which side the boundary is in the element + !! read from database file + !! side_wd = 1--8: only one vertex is on the boundary + !! side_wd = 9--20: only one edge is on the boundary + !! side_wd = 21--26: one face is on the boundary + integer, dimension(:), allocatable :: side_wd + end module meshfem_par diff --git a/src/meshfem3D/rules.mk b/src/meshfem3D/rules.mk index d7e27148b..24cfdaa80 100644 --- a/src/meshfem3D/rules.mk +++ b/src/meshfem3D/rules.mk @@ -59,6 +59,7 @@ meshfem3D_OBJECTS = \ $O/get_MPI_cutplanes_xi.mesh.o \ $O/meshfem3D.mesh.o \ $O/meshfem3D_par.mesh_module.o \ + $O/get_wavefield_discontinuity.mesh.o \ $O/read_mesh_parameter_file.mesh.o \ $O/read_value_mesh_parameters.mesh.o \ $O/save_databases.mesh.o \ diff --git a/src/meshfem3D/save_databases.F90 b/src/meshfem3D/save_databases.F90 index de31b2f64..7562e8986 100644 --- a/src/meshfem3D/save_databases.F90 +++ b/src/meshfem3D/save_databases.F90 @@ -46,6 +46,9 @@ subroutine save_databases(nspec,nglob, & nspec_CPML,is_CPML,CPML_to_spec,CPML_regions, & SAVE_MESH_AS_CUBIT + !! setting up wavefield discontinuity interface + use shared_parameters, only: IS_WAVEFIELD_DISCONTINUITY + implicit none ! number of spectral elements in each block @@ -257,6 +260,7 @@ subroutine save_databases(nspec,nglob, & write(IIN_database) ispec,material_index(1,ispec),material_index(2,ispec),(loc_node(ia),ia = 1,NGNOD) enddo + ! Boundaries ! ! note: check with routine write_boundaries_database() to produce identical output @@ -596,6 +600,10 @@ subroutine save_databases(nspec,nglob, & deallocate(material_index) + ! setting up wavefield discontinuity interface + if (IS_WAVEFIELD_DISCONTINUITY) & + call write_wavefield_discontinuity_database() + ! user output if (myrank == 0) then write(IMAIN,*) ' done mesh files' @@ -623,6 +631,9 @@ subroutine save_output_mesh_files_as_cubit(nspec,nglob, & NMATERIALS,material_properties, & nspec_CPML,CPML_to_spec,CPML_regions + !! setting up wavefield discontinuity interface + use shared_parameters, only: IS_WAVEFIELD_DISCONTINUITY + implicit none integer, parameter :: IIN_database = IIN_DB @@ -778,6 +789,10 @@ subroutine save_output_mesh_files_as_cubit(nspec,nglob, & enddo close(IIN_database) + !! setting up wavefield discontinuity interface + if (IS_WAVEFIELD_DISCONTINUITY) & + call write_wavefield_discontinuity_file() + open(IIN_database,file='MESH/absorbing_surface_file_xmin') write(IIN_database,*) nspec2D_xmin do i = 1,nspec2D_xmin diff --git a/src/shared/count_number_of_sources.f90 b/src/shared/count_number_of_sources.f90 index aba14d4ba..cfde30284 100644 --- a/src/shared/count_number_of_sources.f90 +++ b/src/shared/count_number_of_sources.f90 @@ -37,6 +37,8 @@ subroutine count_number_of_sources(NSOURCES,sources_filename) use shared_parameters, only: USE_FORCE_POINT_SOURCE,USE_EXTERNAL_SOURCE_FILE, & HAS_FINITE_FAULT_SOURCE,NUMBER_OF_SIMULTANEOUS_RUNS + use shared_parameters, only: IS_WAVEFIELD_DISCONTINUITY + implicit none integer,intent(out) :: NSOURCES @@ -125,9 +127,13 @@ subroutine count_number_of_sources(NSOURCES,sources_filename) ! checks if any if (NSOURCES < 1) then - print *,'Error: ',trim(sources_filename),' has ',icounter,'lines, but need ',nlines_per_source, & + ! if there is prescribed wavefield discontinuity, then + ! it is OK to have an empty source file + if (.not. IS_WAVEFIELD_DISCONTINUITY) then + print *,'Error: ',trim(sources_filename),' has ',icounter,'lines, but need ',nlines_per_source, & 'per source... ',NSOURCES - stop 'Error need at least one source in CMTSOLUTION or FORCESOLUTION file' + stop 'Error need at least one source in CMTSOLUTION or FORCESOLUTION file' + endif endif end subroutine count_number_of_sources diff --git a/src/shared/read_parameter_file.F90 b/src/shared/read_parameter_file.F90 index fd56fb6b2..18aa68bfa 100644 --- a/src/shared/read_parameter_file.F90 +++ b/src/shared/read_parameter_file.F90 @@ -711,6 +711,14 @@ subroutine read_parameter_file(BROADCAST_AFTER_READ) !------------------------------------------------------- + ! prescribed wavefield discontinuity on an interface + ! if these parameters do not exist in Par_file, then wavefield + ! is not switched on by default + call read_value_logical(IS_WAVEFIELD_DISCONTINUITY, 'IS_WAVEFIELD_DISCONTINUITY', ier); ier = 0 + if (IS_WAVEFIELD_DISCONTINUITY) write(*,'(a)') 'wavefield discontinuity enabled' + + !------------------------------------------------------- + ! for simultaneous runs from the same batch job call read_value_integer(NUMBER_OF_SIMULTANEOUS_RUNS, 'NUMBER_OF_SIMULTANEOUS_RUNS', ier) if (ier /= 0) then @@ -1057,6 +1065,13 @@ subroutine check_simulation_parameters() MOVIE_VOLUME_STRESS = .false. endif + if (IS_WAVEFIELD_DISCONTINUITY) then + if (GPU_MODE) & + stop 'wavefield discontinuity problem cannot be solved on GPU yet' + if (ATTENUATION) & + stop 'wavefield discontinuity problem cannot be solved with attenuation' + endif + end subroutine check_simulation_parameters ! @@ -1419,6 +1434,9 @@ subroutine broadcast_computed_parameters() call bcast_all_string(FKMODEL_FILE) call bcast_all_singlel(RECIPROCITY_AND_KH_INTEGRAL) + ! wavefield discontinuity + call bcast_all_singlel(IS_WAVEFIELD_DISCONTINUITY) + ! simultaneous runs call bcast_all_singlei(NUMBER_OF_SIMULTANEOUS_RUNS) call bcast_all_singlel(BROADCAST_SAME_MESH_AND_MODEL) diff --git a/src/shared/shared_par.F90 b/src/shared/shared_par.F90 index 36ab72cde..92f50d866 100644 --- a/src/shared/shared_par.F90 +++ b/src/shared/shared_par.F90 @@ -204,6 +204,9 @@ module shared_input_parameters logical :: MESH_A_CHUNK_OF_THE_EARTH logical :: RECIPROCITY_AND_KH_INTEGRAL + ! prescribed wavefield discontinuity on an interface + logical :: IS_WAVEFIELD_DISCONTINUITY = .false. ! if .true. then wavefield discontinuity is turned on (default is false) + end module shared_input_parameters ! diff --git a/src/specfem3D/compute_forces_viscoelastic.F90 b/src/specfem3D/compute_forces_viscoelastic.F90 index 3a5fb47e9..64bec1141 100644 --- a/src/specfem3D/compute_forces_viscoelastic.F90 +++ b/src/specfem3D/compute_forces_viscoelastic.F90 @@ -82,6 +82,11 @@ subroutine compute_forces_viscoelastic(iphase,deltat, & ! PML use pml_par, only: is_CPML,NSPEC_CPML + !! solving wavefield problem with non-split-node scheme + use shared_parameters, only: IS_WAVEFIELD_DISCONTINUITY + use wavefield_discontinuity_solver, only: & + add_displacement_discontinuity_element + ! LTS use specfem_par_lts, only: lts_type_compute_pelem,current_lts_elem,current_lts_boundary_elem @@ -195,6 +200,7 @@ subroutine compute_forces_viscoelastic(iphase,deltat, & !$OMP irregular_element_number,jacobian_regular,xix_regular, & !$OMP displ,veloc,accel, & !$OMP is_CPML,backward_simulation, & +!$OMP IS_WAVEFIELD_DISCONTINUITY, & !$OMP xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,jacobianstore, & !$OMP kappastore,mustore, & !$OMP Kelvin_Voigt_eta,USE_KELVIN_VOIGT_DAMPING, & @@ -312,6 +318,15 @@ subroutine compute_forces_viscoelastic(iphase,deltat, & enddo enddo enddo + !! solving wavefield discontinuity problem with non-split-node scheme + !! add back displacement discontinuity for inner elements of + !! the discontinuity interface + !! note that this is called in adjoint simulation + if (IS_WAVEFIELD_DISCONTINUITY .and. & + ((SIMULATION_TYPE == 1) .or. backward_simulation)) then + call add_displacement_discontinuity_element(ispec, dummyx_loc, & + dummyy_loc, dummyz_loc) + endif endif !------------------------------------------------------------------------------ diff --git a/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 b/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 index 6e770206e..8d3627c0e 100644 --- a/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 +++ b/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 @@ -40,6 +40,10 @@ subroutine compute_forces_viscoelastic_calling() use constants, only: FAULT_SYNCHRONIZE_DISPL_VELOC,FAULT_SYNCHRONIZE_ACCEL use fault_solver_dynamic, only: bc_dynflt_set3d_all,SIMULATION_TYPE_DYN,fault_output_synchronize_GPU,NT_RECORD_LENGTH use fault_solver_kinematic, only: bc_kinflt_set_all,SIMULATION_TYPE_KIN + + !! solving wavefield discontinuity problem with non-split-node scheme + use wavefield_discontinuity_solver, only: & + add_traction_discontinuity, read_wavefield_discontinuity_file implicit none @@ -103,6 +107,13 @@ subroutine compute_forces_viscoelastic_calling() endif endif + !! solving wavefield discontinuity problem with + !! non-split-node scheme + !! note that this is not called in case of adjoint simulation + if (IS_WAVEFIELD_DISCONTINUITY .and. (SIMULATION_TYPE == 1)) then + call read_wavefield_discontinuity_file() + endif + ! distinguishes two runs: for elements in contact with MPI interfaces, and elements within the partitions do iphase = 1,2 @@ -266,6 +277,15 @@ subroutine compute_forces_viscoelastic_calling() ! on GPU call compute_add_sources_viscoelastic_GPU() endif + + !! solving wavefield discontinuity problem with + !! non-split-node scheme + !! traction discontinuity term must be added in iphase = 1 + !! because these terms need to be used in MPI calls + !! note that this is not called in case of adjoint simulation + if (IS_WAVEFIELD_DISCONTINUITY .and. (SIMULATION_TYPE == 1)) then + call add_traction_discontinuity(accel, NGLOB_AB) + endif endif ! iphase ! assemble all the contributions between slices using MPI diff --git a/src/specfem3D/initialize_simulation.F90 b/src/specfem3D/initialize_simulation.F90 index c3baa8bc2..c9617b9d2 100644 --- a/src/specfem3D/initialize_simulation.F90 +++ b/src/specfem3D/initialize_simulation.F90 @@ -345,8 +345,12 @@ subroutine initialize_simulation_check() endif ! check that we have at least one source - if (NSOURCES < 1 .and. (.not. HAS_FINITE_FAULT_SOURCE .and. .not. INVERSE_FWI_FULL_PROBLEM)) & - call exit_MPI(myrank,'need at least one source') + if (NSOURCES < 1 .and. (.not. HAS_FINITE_FAULT_SOURCE .and. .not. INVERSE_FWI_FULL_PROBLEM)) then + ! skip source check in case of wavefield discontinuity problem + if (.not. IS_WAVEFIELD_DISCONTINUITY) then + call exit_MPI(myrank,'need at least one source') + endif + endif ! check simulation type if (SIMULATION_TYPE /= 1 .and. SIMULATION_TYPE /= 2 .and. SIMULATION_TYPE /= 3) & diff --git a/src/specfem3D/iterate_time.F90 b/src/specfem3D/iterate_time.F90 index 5f16e06c0..5c413af04 100644 --- a/src/specfem3D/iterate_time.F90 +++ b/src/specfem3D/iterate_time.F90 @@ -39,6 +39,12 @@ subroutine iterate_time() ! hdf5 i/o server use io_server_hdf5, only: do_io_start_idle,pass_info_to_io + !! solving wavefield discontinuity problem with non-split-node scheme + use wavefield_discontinuity_solver, only: & + read_mesh_databases_wavefield_discontinuity, & + open_wavefield_discontinuity_file, & + finalize_wavefield_discontinuity + implicit none ! for EXACT_UNDOING_TO_DISK @@ -207,6 +213,13 @@ subroutine iterate_time() ! get MPI starting time_start = wtime() + !! solving wavefield discontinuity problem with + !! non-split-node scheme + if (IS_WAVEFIELD_DISCONTINUITY) then + call read_mesh_databases_wavefield_discontinuity() + call open_wavefield_discontinuity_file() + endif + ! LTS if (LTS_MODE) then ! LTS steps through its own time iterations - for now @@ -329,6 +342,12 @@ subroutine iterate_time() ! enddo ! end of main time loop + !! solving wavefield discontinuity problem with + !! non-split-node scheme + if (IS_WAVEFIELD_DISCONTINUITY) then + call finalize_wavefield_discontinuity() + endif + ! goto point for LTS to finish time loop 777 continue diff --git a/src/specfem3D/rules.mk b/src/specfem3D/rules.mk index 76595154a..d1d81a9df 100644 --- a/src/specfem3D/rules.mk +++ b/src/specfem3D/rules.mk @@ -45,6 +45,7 @@ specfem3D_TARGETS = \ specfem3D_OBJECTS = \ $O/specfem3D_par.spec_module.o \ $O/asdf_data.spec_module.o \ + $O/wavefield_discontinuity_solver.spec.o \ $O/assemble_MPI_vector.spec.o \ $O/check_stability.spec.o \ $O/comp_source_time_function.spec.o \ @@ -199,6 +200,7 @@ specfem3D_MODULES = \ $(FC_MODDIR)/specfem_par_noise.$(FC_MODEXT) \ $(FC_MODDIR)/specfem_par_lts.$(FC_MODEXT) \ $(FC_MODDIR)/user_noise_distribution.$(FC_MODEXT) \ + $(FC_MODDIR)/wavefield_discontinuity_solver.$(FC_MODEXT) \ $(EMPTY_MACRO) @@ -358,6 +360,11 @@ $O/initialize_simulation.spec.o: $O/adios_manager.shared_adios_module.o ## ASDF compilation $O/write_output_ASDF.spec.o: $O/asdf_data.spec_module.o +## wavefield discontinuity +$O/compute_forces_viscoelastic_calling_routine.spec.o: $O/wavefield_discontinuity_solver.spec.o +$O/compute_forces_viscoelastic.spec.o: $O/wavefield_discontinuity_solver.spec.o +$O/iterate_time.spec.o: $O/wavefield_discontinuity_solver.spec.o + ## kdtree $O/locate_point.spec.o: $O/search_kdtree.shared.o $O/setup_sources_receivers.spec.o: $O/search_kdtree.shared.o diff --git a/src/specfem3D/wavefield_discontinuity_solver.f90 b/src/specfem3D/wavefield_discontinuity_solver.f90 new file mode 100644 index 000000000..e13886ffa --- /dev/null +++ b/src/specfem3D/wavefield_discontinuity_solver.f90 @@ -0,0 +1,181 @@ +!! Solving the wavefield discontinuity problem with a non-split-node +!! scheme +!! Tianshi Liu, 2023.5 +module wavefield_discontinuity_solver + use constants, only: CUSTOM_REAL + + !! ispec_to_elem_wd(NSPEC_AB) + !! ispec_to_elem_wd(ispec) = ispec_wd (0 if element not belong to boundary) + !! read from solver database + integer, dimension(:), allocatable :: ispec_to_elem_wd + + !! number of distinct gll points on the boundary + !! read from solver database + integer :: nglob_wd + + !! number of elements on the inner side of the boundary + !! read from solver database + integer :: nspec_wd + + !! ibool_wd(NGLLX, NGLLY, NGLLZ, nspec_wd) + !! ibool_wd(i,j,k,ispec_wd) = iglob_wd (0 if point not on boundary) + !! read from solver database + integer, dimension(:,:,:,:), allocatable :: ibool_wd + + !! boundary_to_iglob_wd(nglob_wd) + !! boundary_to_iglob_wd(iglob_wd) = iglob + !! read from solver database + integer, dimension(:), allocatable :: boundary_to_iglob_wd + + !! mass_in_wd(nglob_wd) + !! mass matrix on the inner side of the boundary + !! note that it is not assembled over processors + !! read from solver database + real(kind=CUSTOM_REAL), dimension(:), allocatable :: mass_in_wd + + !! number of faces on the boundary + !! read from solver database + integer :: nfaces_wd + + !! face_ijk_wd(NDIM, NGLLSQUARE, nfaces_wd) + !! read from solver database + integer, dimension(:,:,:), allocatable :: face_ijk_wd + + !! face_ispec_wd(nfaces_wd) + !! read from solver database + integer, dimension(:), allocatable :: face_ispec_wd + + !! face_normal_wd(NDIM, NGLLSQUARE, nfaces_wd) + !! read from solver database + real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: face_normal_wd + + !! face_jacobian2Dw_wd(NGLLSQUARE, nfaces_wd) + !! read from solver database + real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: face_jacobian2dw_wd + + !! displ_wd(NDIM, nglob_wd) + !! displacement discontinuity condition at current time step + real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: displ_wd + + !! accel_wd(NDIM, nglob_wd) + !! acceleration discontinuity condition at current time step + real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_wd + + !! traction_wd(NDIM, NGLLSQUARE, nfaces_wd) + !! traction discontinuity condition at current time step + real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: traction_wd + +contains + + subroutine read_mesh_databases_wavefield_discontinuity() + use constants, only: IFILE_WAVEFIELD_DISCONTINUITY, & + FNAME_WAVEFIELD_DISCONTINUITY_DATABASE + use specfem_par, only: CUSTOM_REAL, prname, & + NSPEC_AB, NGLLX, NGLLY, NGLLZ, NDIM, NGLLSQUARE + implicit none + open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & + file=trim(prname)//trim(FNAME_WAVEFIELD_DISCONTINUITY_DATABASE), & + action='read', form='unformatted') + + allocate(ispec_to_elem_wd(NSPEC_AB)) + read(IFILE_WAVEFIELD_DISCONTINUITY) ispec_to_elem_wd + read(IFILE_WAVEFIELD_DISCONTINUITY) nglob_wd + read(IFILE_WAVEFIELD_DISCONTINUITY) nspec_wd + allocate(ibool_wd(NGLLX, NGLLY, NGLLZ, nspec_wd), & + boundary_to_iglob_wd(nglob_wd), & + mass_in_wd(nglob_wd)) + read(IFILE_WAVEFIELD_DISCONTINUITY) ibool_wd + read(IFILE_WAVEFIELD_DISCONTINUITY) boundary_to_iglob_wd + read(IFILE_WAVEFIELD_DISCONTINUITY) mass_in_wd + read(IFILE_WAVEFIELD_DISCONTINUITY) nfaces_wd + allocate(face_ijk_wd(NDIM, NGLLSQUARE, nfaces_wd), & + face_ispec_wd(nfaces_wd), & + face_normal_wd(NDIM, NGLLSQUARE, nfaces_wd), & + face_jacobian2Dw_wd(NGLLSQUARE, nfaces_wd)) + read(IFILE_WAVEFIELD_DISCONTINUITY) face_ijk_wd + read(IFILE_WAVEFIELD_DISCONTINUITY) face_ispec_wd + read(IFILE_WAVEFIELD_DISCONTINUITY) face_normal_wd + read(IFILE_WAVEFIELD_DISCONTINUITY) face_jacobian2Dw_wd + allocate(displ_wd(NDIM, nglob_wd), & + accel_wd(NDIM, nglob_wd), & + traction_wd(NDIM, NGLLSQUARE, nfaces_wd)) + end subroutine read_mesh_databases_wavefield_discontinuity + + subroutine open_wavefield_discontinuity_file() + use specfem_par, only: prname + use constants, only: IFILE_WAVEFIELD_DISCONTINUITY + implicit none + open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & + file=trim(prname)//'wavefield_discontinuity.bin', & + status='old',action='read',form='unformatted') + end subroutine open_wavefield_discontinuity_file + + subroutine read_wavefield_discontinuity_file() + use constants, only: IFILE_WAVEFIELD_DISCONTINUITY + implicit none + read(IFILE_WAVEFIELD_DISCONTINUITY) displ_wd + read(IFILE_WAVEFIELD_DISCONTINUITY) accel_wd + read(IFILE_WAVEFIELD_DISCONTINUITY) traction_wd + end subroutine read_wavefield_discontinuity_file + + subroutine finalize_wavefield_discontinuity() + use constants, only: IFILE_WAVEFIELD_DISCONTINUITY + implicit none + close(IFILE_WAVEFIELD_DISCONTINUITY) + deallocate(ispec_to_elem_wd, ibool_wd, boundary_to_iglob_wd, mass_in_wd, & + face_ijk_wd, face_ispec_wd, face_normal_wd, face_jacobian2Dw_wd, & + displ_wd, accel_wd, traction_wd) + end subroutine finalize_wavefield_discontinuity + + subroutine add_displacement_discontinuity_element(ispec, dummyx_loc, & + dummyy_loc, dummyz_loc) + use specfem_par, only: CUSTOM_REAL,NGLLX, NGLLY, NGLLZ + implicit none + integer, intent(in) :: ispec + real(kind=CUSTOM_REAL), dimension(NGLLX, NGLLY, NGLLZ) :: dummyx_loc, & + dummyy_loc, dummyz_loc + integer :: ispec_wd, i, j, k, iglob_wd + ispec_wd = ispec_to_elem_wd(ispec) + if (ispec_wd /= 0) then + do k=1,NGLLZ + do j=1,NGLLY + do i=1,NGLLX + iglob_wd = ibool_wd(i,j,k,ispec_wd) + if (iglob_wd /= 0) then + dummyx_loc(i,j,k) = dummyx_loc(i,j,k) + displ_wd(1, iglob_wd) + dummyy_loc(i,j,k) = dummyy_loc(i,j,k) + displ_wd(2, iglob_wd) + dummyz_loc(i,j,k) = dummyz_loc(i,j,k) + displ_wd(3, iglob_wd) + endif + enddo + enddo + enddo + endif + end subroutine add_displacement_discontinuity_element + + subroutine add_traction_discontinuity(accel, nglob) + use specfem_par, only: CUSTOM_REAL, NGLLX, NGLLY, NGLLZ, NGLLSQUARE, & + ibool, NDIM + !use specfem_par_elastic, only: accel + implicit none + integer :: iglob_wd, iglob, ispec, i, j, k, iface_wd, igll, nglob + real(kind=CUSTOM_REAL) :: jacobianw + real(kind=CUSTOM_REAL) :: accel(NDIM, nglob) + do iglob_wd = 1, nglob_wd + iglob = boundary_to_iglob_wd(iglob_wd) + accel(:,iglob) = accel(:,iglob) - & + accel_wd(:,iglob_wd) * mass_in_wd(iglob_wd) + enddo + do iface_wd = 1, nfaces_wd + do igll = 1, NGLLSQUARE + i = face_ijk_wd(1, igll, iface_wd) + j = face_ijk_wd(2, igll, iface_wd) + k = face_ijk_wd(3, igll, iface_wd) + ispec = face_ispec_wd(iface_wd) + iglob = ibool(i,j,k,ispec) + jacobianw = face_jacobian2Dw_wd(igll, iface_wd) + accel(:,iglob) = accel(:,iglob) + & + traction_wd(:,igll,iface_wd) * jacobianw + enddo + enddo + end subroutine add_traction_discontinuity +end module wavefield_discontinuity_solver