[cig-commits] r11867 - seismo/3D/SPECFEM3D_GLOBE/trunk
liuqy at geodynamics.org
liuqy at geodynamics.org
Tue Apr 29 11:02:54 PDT 2008
Author: liuqy
Date: 2008-04-29 11:02:53 -0700 (Tue, 29 Apr 2008)
New Revision: 11867
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_arrays_source.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
Log:
Add source frechet derivatives output for adjoint simulations
(SIMULATION_TYPE == 2). Output file is in OUTPUT_FILES/directory
with names like src_frechet.00001
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_arrays_source.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_arrays_source.f90 2008-04-25 20:44:32 UTC (rev 11866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_arrays_source.f90 2008-04-29 18:02:53 UTC (rev 11867)
@@ -328,4 +328,185 @@
end subroutine comp_subarrays_adjoint_src
+! =======================================================================
+! compute the integrated derivatives of source parameters (M_jk and X_s)
+
+subroutine compute_adj_source_frechet(displ_s,Mxx,Myy,Mzz,Mxy,Mxz,Myz,eps_s,eps_m_s, &
+ hxir,hetar,hgammar,hpxir,hpetar,hpgammar, hprime_xx,hprime_yy,hprime_zz, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+
+ implicit none
+
+ include 'constants.h'
+
+ ! input
+ real(kind=CUSTOM_REAL) :: displ_s(NDIM,NGLLX,NGLLY,NGLLZ)
+ double precision :: Mxx, Myy, Mzz, Mxy, Mxz, Myz
+ ! output
+ real(kind=CUSTOM_REAL) :: eps_s(NDIM,NDIM), eps_m_s(NDIM)
+
+ ! auxilliary
+ double precision :: hxir(NGLLX), hetar(NGLLY), hgammar(NGLLZ), &
+ hpxir(NGLLX),hpetar(NGLLY),hpgammar(NGLLZ)
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+! local variables
+ real(kind=CUSTOM_REAL) :: tempx1l,tempx2l,tempx3l, tempy1l,tempy2l,tempy3l, &
+ tempz1l,tempz2l,tempz3l, hp1, hp2, hp3, &
+ xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl, &
+ duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl, &
+ xix_s,xiy_s,xiz_s,etax_s,etay_s,etaz_s,gammax_s,gammay_s,gammaz_s, &
+ hlagrange_xi, hlagrange_eta, hlagrange_gamma, hlagrange
+
+ real(kind=CUSTOM_REAL) :: eps(NDIM,NDIM), eps_array(NDIM,NDIM,NGLLX,NGLLY,NGLLZ), &
+ eps_m_array(NGLLX,NGLLY,NGLLZ)
+
+ integer i,j,k,l
+
+
+! first compute the strain at all the GLL points of the source element
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ tempx1l = 0._CUSTOM_REAL
+ tempx2l = 0._CUSTOM_REAL
+ tempx3l = 0._CUSTOM_REAL
+
+ tempy1l = 0._CUSTOM_REAL
+ tempy2l = 0._CUSTOM_REAL
+ tempy3l = 0._CUSTOM_REAL
+
+ tempz1l = 0._CUSTOM_REAL
+ tempz2l = 0._CUSTOM_REAL
+ tempz3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ tempx1l = tempx1l + displ_s(1,l,j,k)*hp1
+ tempy1l = tempy1l + displ_s(2,l,j,k)*hp1
+ tempz1l = tempz1l + displ_s(3,l,j,k)*hp1
+
+ hp2 = hprime_yy(j,l)
+ tempx2l = tempx2l + displ_s(1,i,l,k)*hp2
+ tempy2l = tempy2l + displ_s(2,i,l,k)*hp2
+ tempz2l = tempz2l + displ_s(3,i,l,k)*hp2
+
+ hp3 = hprime_zz(k,l)
+ tempx3l = tempx3l + displ_s(1,i,j,l)*hp3
+ tempy3l = tempy3l + displ_s(2,i,j,l)*hp3
+ tempz3l = tempz3l + displ_s(3,i,j,l)*hp3
+ enddo
+
+! dudx
+ xixl = xix(i,j,k)
+ xiyl = xiy(i,j,k)
+ xizl = xiz(i,j,k)
+ etaxl = etax(i,j,k)
+ etayl = etay(i,j,k)
+ etazl = etaz(i,j,k)
+ gammaxl = gammax(i,j,k)
+ gammayl = gammay(i,j,k)
+ gammazl = gammaz(i,j,k)
+
+ duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+! strain eps_jk
+ eps(1,1) = duxdxl
+ eps(1,2) = (duxdyl + duydxl) / 2
+ eps(1,3) = (duxdzl + duzdxl) / 2
+ eps(2,2) = duydyl
+ eps(2,3) = (duydzl + duzdyl) / 2
+ eps(3,3) = duzdzl
+ eps(2,1) = eps(1,2)
+ eps(3,1) = eps(1,3)
+ eps(3,2) = eps(2,3)
+
+ eps_array(:,:,i,j,k) = eps(:,:)
+
+! Mjk eps_jk
+ eps_m_array(i,j,k) = Mxx * eps(1,1) + Myy * eps(2,2) + Mzz * eps(3,3) + &
+ 2 * (Mxy * eps(1,2) + Mxz * eps(1,3) + Myz * eps(2,3))
+
+ enddo
+ enddo
+ enddo
+
+ ! interpolate the strain eps_s(:,:) from eps_array(:,:,i,j,k)
+ eps_s = 0.
+ xix_s = 0.; xiy_s = 0.; xiz_s = 0.
+ etax_s = 0.; etay_s = 0.; etaz_s = 0.
+ gammax_s = 0.; gammay_s = 0.; gammaz_s = 0.
+
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ hlagrange = hxir(i)*hetar(j)*hgammar(k)
+
+ eps_s(1,1) = eps_s(1,1) + eps_array(1,1,i,j,k)*hlagrange
+ eps_s(1,2) = eps_s(1,2) + eps_array(1,2,i,j,k)*hlagrange
+ eps_s(1,3) = eps_s(1,3) + eps_array(1,3,i,j,k)*hlagrange
+ eps_s(2,2) = eps_s(2,2) + eps_array(2,2,i,j,k)*hlagrange
+ eps_s(2,3) = eps_s(2,3) + eps_array(2,3,i,j,k)*hlagrange
+ eps_s(3,3) = eps_s(3,3) + eps_array(3,3,i,j,k)*hlagrange
+
+ xix_s = xix_s + xix(i,j,k)*hlagrange
+ xiy_s = xiy_s + xiy(i,j,k)*hlagrange
+ xiz_s = xiz_s + xiz(i,j,k)*hlagrange
+ etax_s = etax_s + etax(i,j,k)*hlagrange
+ etay_s = etay_s + etay(i,j,k)*hlagrange
+ etaz_s = etaz_s + etaz(i,j,k)*hlagrange
+ gammax_s = gammax_s + gammax(i,j,k)*hlagrange
+ gammay_s = gammay_s + gammay(i,j,k)*hlagrange
+ gammaz_s = gammaz_s + gammaz(i,j,k)*hlagrange
+
+ enddo
+ enddo
+ enddo
+
+! for completion purpose, not used in specfem3D.f90
+ eps_s(2,1) = eps_s(1,2)
+ eps_s(3,1) = eps_s(1,3)
+ eps_s(3,2) = eps_s(2,3)
+
+! compute the gradient of M_jk * eps_jk, and then interpolate it
+
+ eps_m_s = 0.
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ hlagrange_xi = hpxir(i)*hetar(j)*hgammar(k)
+ hlagrange_eta = hxir(i)*hpetar(j)*hgammar(k)
+ hlagrange_gamma = hxir(i)*hetar(j)*hpgammar(k)
+
+ eps_m_s(1) = eps_m_s(1) + eps_m_array(i,j,k) * (hlagrange_xi * xix_s &
+ + hlagrange_eta * etax_s + hlagrange_gamma * gammax_s)
+ eps_m_s(2) = eps_m_s(2) + eps_m_array(i,j,k) * (hlagrange_xi * xiy_s &
+ + hlagrange_eta * etay_s + hlagrange_gamma * gammay_s)
+ eps_m_s(3) = eps_m_s(3) + eps_m_array(i,j,k) * (hlagrange_xi * xiz_s &
+ + hlagrange_eta * etaz_s + hlagrange_gamma * gammaz_s)
+
+ enddo
+ enddo
+ enddo
+
+end subroutine compute_adj_source_frechet
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2008-04-25 20:44:32 UTC (rev 11866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2008-04-29 18:02:53 UTC (rev 11867)
@@ -499,7 +499,7 @@
! Newmark time scheme parameters and non-dimensionalization
real(kind=CUSTOM_REAL) time,deltat,deltatover2,deltatsqover2
- double precision scale_t,scale_displ,scale_veloc
+ double precision scale_t,scale_displ,scale_veloc,scale_mass
! ADJOINT
real(kind=CUSTOM_REAL) b_additional_term,b_force_normal_comp
@@ -588,7 +588,13 @@
integer NSTEP_SUB_ADJ,it_sub_adj,iadj_block ! to read input in chunks
integer, dimension(:,:), allocatable :: iadjsrc ! to read input in chunks
integer, dimension(:), allocatable :: iadjsrc_len,iadj_vec
+! source frechet derivatives
+ real(kind=CUSTOM_REAL) :: displ_s(NDIM,NGLLX,NGLLY,NGLLZ), eps_s(NDIM,NDIM), eps_m_s(NDIM), stf_deltat
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: moment_der
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: sloc_der
+ double precision, dimension(:,:), allocatable :: hpxir_store,hpetar_store,hpgammar_store
+
! seismograms
integer it_begin,it_end,nit_written
double precision uxd, uyd, uzd, eps_trace,dxx,dyy,dxy,dxz,dyz,eps_loc(NDIM,NDIM), eps_loc_new(NDIM,NDIM)
@@ -1988,6 +1994,9 @@
hgammar_store(irec_local,:) = hgammar(:)
enddo
else
+ allocate(hpxir_store(nrec_local,NGLLX))
+ allocate(hpetar_store(nrec_local,NGLLY))
+ allocate(hpgammar_store(nrec_local,NGLLZ))
do irec_local = 1,nrec_local
irec = number_receiver_global(irec_local)
call lagrange_any(xi_source(irec),NGLLX,xigll,hxir,hpxir)
@@ -1996,6 +2005,9 @@
hxir_store(irec_local,:) = hxir(:)
hetar_store(irec_local,:) = hetar(:)
hgammar_store(irec_local,:) = hgammar(:)
+ hpxir_store(irec_local,:) = hpxir(:)
+ hpetar_store(irec_local,:) = hpetar(:)
+ hpgammar_store(irec_local,:) = hpgammar(:)
enddo
endif
@@ -2439,8 +2451,12 @@
allocate(seismograms(NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
if(ier /= 0) stop 'error while allocating seismograms'
else
- allocate(seismograms(9,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
+ allocate(seismograms(NDIM*NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
if(ier /= 0) stop 'error while allocating seismograms'
+ ! allocate Frechet derivatives array
+ allocate(moment_der(NDIM,NDIM,nrec_local),sloc_der(NDIM,nrec_local))
+ moment_der = 0._CUSTOM_REAL
+ sloc_der = 0._CUSTOM_REAL
endif
! initialize seismograms
seismograms(:,:,:) = 0._CUSTOM_REAL
@@ -4375,6 +4391,8 @@
dxy = dxy + dble(epsilondev_crust_mantle(3,i,j,k,ispec_selected_source(irec)))*hlagrange
dxz = dxz + dble(epsilondev_crust_mantle(4,i,j,k,ispec_selected_source(irec)))*hlagrange
dyz = dyz + dble(epsilondev_crust_mantle(5,i,j,k,ispec_selected_source(irec)))*hlagrange
+
+ displ_s(:,i,j,k) = displ_crust_mantle(:,iglob)
enddo
enddo
@@ -4391,7 +4409,7 @@
eps_loc(3,2) = dyz
eps_loc_new(:,:) = eps_loc(:,:)
-! rotate to the local cartesian coordinates (e-n-z)
+! rotate to the local cartesian coordinates (e-n-z): eps_new=P*eps*P'
eps_loc_new(:,:) = matmul(matmul(nu_source(:,:,irec),eps_loc(:,:)), transpose(nu_source(:,:,irec)))
! distinguish between single and double precision for reals
@@ -4415,6 +4433,22 @@
nu_source(:,2,irec)*uyd + nu_source(:,3,irec)*uzd)
endif
+! frechet derviatives of the source
+ ispec = ispec_selected_source(irec)
+
+ call compute_adj_source_frechet(displ_s,Mxx(irec),Myy(irec),Mzz(irec),Mxy(irec),Mxz(irec),Myz(irec),eps_s,eps_m_s, &
+ hxir_store(irec_local,:),hetar_store(irec_local,:),hgammar_store(irec_local,:), &
+ hpxir_store(irec_local,:),hpetar_store(irec_local,:),hpgammar_store(irec_local,:),hprime_xx,hprime_yy,hprime_zz, &
+ xix_crust_mantle(:,:,:,ispec),xiy_crust_mantle(:,:,:,ispec),xiz_crust_mantle(:,:,:,ispec), &
+ etax_crust_mantle(:,:,:,ispec),etay_crust_mantle(:,:,:,ispec),etaz_crust_mantle(:,:,:,ispec), &
+ gammax_crust_mantle(:,:,:,ispec),gammay_crust_mantle(:,:,:,ispec),gammaz_crust_mantle(:,:,:,ispec))
+
+ stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-t_cmt(irec),hdur_gaussian(irec))
+ stf_deltat = stf * deltat
+ moment_der(:,:,irec_local) = moment_der(:,:,irec_local) + eps_s(:,:) * stf_deltat
+ sloc_der(:,irec_local) = sloc_der(:,irec_local) + eps_m_s(:) * stf_deltat
+
+
else if (SIMULATION_TYPE == 3) then
do k = 1,NGLLZ
@@ -5110,6 +5144,43 @@
endif
+! save source derivatives for adjoint simulations
+ if (SIMULATION_TYPE == 2 .and. nrec_local > 0) then
+ scale_mass = RHOAV * (R_EARTH**3)
+
+ do irec_local = 1, nrec_local
+! rotate and scale the location derivatives to correspond to dn,de,dz
+ sloc_der(:,irec_local) = matmul(nu_source(:,:,irec),sloc_der(:,irec_local)) * scale_displ * scale_t
+! rotate scale the moment derivatives to correspond to M[n,e,z][n,e,z]
+ moment_der(:,:,irec_local) = matmul(matmul(nu_source(:,:,irec),moment_der(:,:,irec_local)),&
+ transpose(nu_source(:,:,irec))) * scale_t ** 3 / scale_mass
+
+ write(outputname,'(a,i5.5)') 'OUTPUT_FILES/src_frechet.',number_receiver_global(irec_local)
+ open(unit=27,file=trim(outputname),status='unknown')
+!
+! r -> z, theta -> -n, phi -> e, plus factor 2 for Mrt,Mrp,Mtp, and 1e-7 to dyne.cm
+! Mrr = Mzz
+! Mtt = Mnn
+! Mpp = Mee
+! Mrt = -Mzn
+! Mrp = Mze
+! Mtp = -Mne
+! minus sign for sloc_der(3,irec_local) to get derivative for depth instead of radius
+
+ write(27,'(g16.5)') moment_der(3,3,irec_local) * 1e-7
+ write(27,'(g16.5)') moment_der(1,1,irec_local) * 1e-7
+ write(27,'(g16.5)') moment_der(2,2,irec_local) * 1e-7
+ write(27,'(g16.5)') -2*moment_der(1,3,irec_local) * 1e-7
+ write(27,'(g16.5)') 2*moment_der(2,3,irec_local) * 1e-7
+ write(27,'(g16.5)') -2*moment_der(1,2,irec_local) * 1e-7
+ write(27,'(g16.5)') sloc_der(2,irec_local)
+ write(27,'(g16.5)') sloc_der(1,irec_local)
+ write(27,'(g16.5)') -sloc_der(3,irec_local)
+ close(27)
+ enddo
+ endif
+
+
! if running on MareNostrum in Barcelona
if(RUN_ON_MARENOSTRUM_BARCELONA) then
More information about the cig-commits
mailing list