[cig-commits] r19844 - in seismo/3D/ADJOINT_TOMO/iterate_adj: UTILS cluster cluster/model_slice cluster/model_vp_vs cluster/model_vp_vs/src cluster/scripts cluster/smooth cluster/smooth_update cluster/sum_kernel cluster/sum_kernel/src matlab
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Wed Mar 21 18:34:38 PDT 2012
Author: danielpeter
Date: 2012-03-21 18:34:38 -0700 (Wed, 21 Mar 2012)
New Revision: 19844
Removed:
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_slice/exit_mpi.o
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_slice/read_basin_topo_bathy_file.o
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_slice/sem_model_slice.o
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_slice/utm_geo.o
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/xcombine_vol_data
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/smooth_sem_fun
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/xcombine_vol_data
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth_update/smooth_sem_fun
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth_update/xcombine_vol_data
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/kernels_run
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/sum_kernels
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/xcombine_vol_data
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/kernel_combine.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/README
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso_cg.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso_iso.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/scripts/README
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/sum_kernels_globe.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/sum_preconditioned_kernels_globe.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/README
Log:
updates kernel summation and global model updates; removes compiled binaries and object files from repository
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/kernel_combine.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/kernel_combine.pl 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/kernel_combine.pl 2012-03-22 01:34:38 UTC (rev 19844)
@@ -9,6 +9,7 @@
# EXAMPLE:
# ~/UTILS/kernel_combine.pl 5 7 mu_all_kernel_set
#
+# combines pdf files for a kernel plot
#==========================================================
if (@ARGV < 3) {die("Usage: kernel_combine.pl imodel_min imodel_max suffix\n")}
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/README 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/README 2012-03-22 01:34:38 UTC (rev 19844)
@@ -5,7 +5,9 @@
/ADJOINT_TOMO/iterate_adj/cluster/README
This is designed to provide guidance for doing an inversion using event kernels
-derived from adjoint methods. Programs are contained for the basic operations needed for a tomographic inversion based on adjoint methods. However, several minor adjustments need to be made for each user.
+derived from adjoint methods. Programs are contained for the basic operations
+needed for a tomographic inversion based on adjoint methods.
+However, several minor adjustments need to be made for each user.
Email me (carltape at fas.harvard.edu) with suggestions for improvements.
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_slice/exit_mpi.o
===================================================================
(Binary files differ)
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_slice/read_basin_topo_bathy_file.o
===================================================================
(Binary files differ)
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_slice/sem_model_slice.o
===================================================================
(Binary files differ)
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_slice/utm_geo.o
===================================================================
(Binary files differ)
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model
===================================================================
(Binary files differ)
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso.f90 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso.f90 2012-03-22 01:34:38 UTC (rev 19844)
@@ -79,7 +79,7 @@
! volume
real(kind=CUSTOM_REAL), dimension(NGLOB) :: x, y, z
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
- integer, dimension(NSPEC) :: idoubling
+ logical, dimension(NSPEC) :: ispec_is_tiso
! gradient vector norm ( v^T * v )
real(kind=CUSTOM_REAL) :: norm_bulk,norm_betav,norm_betah,norm_eta
@@ -159,8 +159,8 @@
alphah1 = 0._CUSTOM_REAL
! do not use transverse isotropy except if element is between d220 and Moho
- if(.not. ( idoubling(ispec)==IFLAG_220_80 .or. idoubling(ispec)==IFLAG_80_MOHO) ) then
-
+! if(.not. ( idoubling(ispec)==IFLAG_220_80 .or. idoubling(ispec)==IFLAG_80_MOHO) ) then
+ if( .not. ispec_is_tiso(ispec) ) then
! isotropic model update
! no eta perturbation, since eta = 1 in isotropic media
@@ -349,6 +349,8 @@
use model_update_tiso
implicit none
+ integer, dimension(:), allocatable :: idummy
+ logical, dimension(:), allocatable :: ldummy
! vpv model
write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vpv.bin'
@@ -440,6 +442,10 @@
print*
endif
+ ! allocates temporary array
+ allocate(idummy(NSPEC),ldummy(NSPEC),stat=ier)
+ if( ier /= 0 ) stop 'error allocating ldummy'
+
! global addressing
write(m_file,'(a,i6.6,a)') 'topo/proc',myrank,'_reg1_solver_data_2.bin'
open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
@@ -451,9 +457,13 @@
read(11) y(1:nglob)
read(11) z(1:nglob)
read(11) ibool(:,:,:,1:nspec)
- read(11) idoubling(1:nspec)
+ read(11) idummy(1:nspec) ! idoubling
+ read(11) ldummy(1:nspec) ! is_on_a_slice_edge
+ read(11) ispec_is_tiso(1:nspec)
close(11)
+ deallocate(idummy,ldummy)
+
end subroutine read_model
!
@@ -946,11 +956,11 @@
integral_betah,integral_eta
real(kind=CUSTOM_REAL) :: volume_glob,volume_glob_sum
! root-mean square values
- real(kind=CUSTOM_REAL) :: rms_vpv,rms_vph,rms_vsv,rms_vsh,rms_eta,rms_rho
+ real(kind=CUSTOM_REAL) :: rms_vpv,rms_vph,rms_vsv,rms_vsh,rms_eta,rms_rho
real(kind=CUSTOM_REAL) :: rms_vpv_sum,rms_vph_sum,rms_vsv_sum,rms_vsh_sum, &
rms_eta_sum,rms_rho_sum
real(kind=CUSTOM_REAL) :: dvpv,dvph,dvsv,dvsh,deta,drho
-
+
! Gauss-Lobatto-Legendre points of integration and weights
double precision, dimension(NGLLX) :: xigll, wxgll
double precision, dimension(NGLLY) :: yigll, wygll
@@ -1032,9 +1042,9 @@
rms_vpv = 0._CUSTOM_REAL
rms_vph = 0._CUSTOM_REAL
rms_vsv = 0._CUSTOM_REAL
- rms_vsh = 0._CUSTOM_REAL
- rms_eta = 0._CUSTOM_REAL
- rms_rho = 0._CUSTOM_REAL
+ rms_vsh = 0._CUSTOM_REAL
+ rms_eta = 0._CUSTOM_REAL
+ rms_rho = 0._CUSTOM_REAL
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -1081,11 +1091,11 @@
call exit_MPI(myrank,'error NaN')
endif
- ! root-mean square
+ ! root-mean square
! integrates relative perturbations ( dv / v using logarithm ) squared
dvpv = log( model_vpv_new(i,j,k,ispec) / model_vpv(i,j,k,ispec) ) ! alphav
rms_vpv = rms_vpv + volumel * dvpv*dvpv
-
+
dvph = log( model_vph_new(i,j,k,ispec) / model_vph(i,j,k,ispec) ) ! alphah
rms_vph = rms_vph + volumel * dvph*dvph
@@ -1097,7 +1107,7 @@
deta = log( model_eta_new(i,j,k,ispec) / model_eta(i,j,k,ispec) ) ! eta
rms_eta = rms_eta + volumel * deta*deta
-
+
drho = log( model_rho_new(i,j,k,ispec) / model_rho(i,j,k,ispec) ) ! rho
rms_rho = rms_rho + volumel * drho*drho
@@ -1144,7 +1154,7 @@
print*
endif
- ! root-mean square
+ ! root-mean square
call mpi_reduce(rms_vpv,rms_vpv_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(rms_vph,rms_vph_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(rms_vsv,rms_vsv_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso_cg.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso_cg.f90 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso_cg.f90 2012-03-22 01:34:38 UTC (rev 19844)
@@ -105,7 +105,7 @@
! volume
real(kind=CUSTOM_REAL), dimension(NGLOB) :: x, y, z
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
- integer, dimension(NSPEC) :: idoubling
+ logical, dimension(NSPEC) :: ispec_is_tiso
! gradient vector norm ( v^T * v )
real(kind=CUSTOM_REAL) :: norm_bulk,norm_betav,norm_betah,norm_eta
@@ -190,8 +190,8 @@
alphah1 = 0._CUSTOM_REAL
! do not use transverse isotropy except if element is between d220 and Moho
- if(.not. ( idoubling(ispec)==IFLAG_220_80 .or. idoubling(ispec)==IFLAG_80_MOHO) ) then
-
+! if(.not. ( idoubling(ispec)==IFLAG_220_80 .or. idoubling(ispec)==IFLAG_80_MOHO) ) then
+ if( .not. ispec_is_tiso(ispec) ) then
! isotropic model update
! no eta perturbation, since eta = 1 in isotropic media
@@ -386,6 +386,8 @@
use model_update_cg
implicit none
+ integer, dimension(:), allocatable :: idummy
+ logical, dimension(:), allocatable :: ldummy
! vpv model
write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vpv.bin'
@@ -477,6 +479,10 @@
print*
endif
+ ! allocates temporary array
+ allocate(idummy(NSPEC),ldummy(NSPEC),stat=ier)
+ if( ier /= 0 ) stop 'error allocating ldummy'
+
! global addressing
write(m_file,'(a,i6.6,a)') 'topo/proc',myrank,'_reg1_solver_data_2.bin'
open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
@@ -488,9 +494,13 @@
read(11) y(1:nglob)
read(11) z(1:nglob)
read(11) ibool(:,:,:,1:nspec)
- read(11) idoubling(1:nspec)
+ read(11) idummy(1:nspec) ! idoubling
+ read(11) ldummy(1:nspec) ! is_on_a_slice_edge
+ read(11) ispec_is_tiso(1:nspec)
close(11)
+ deallocate(idummy,ldummy)
+
end subroutine read_model
!
@@ -783,6 +793,7 @@
real(kind=CUSTOM_REAL) :: minmax(4),depthmax(2),depthmax_radius(2),max
real(kind=CUSTOM_REAL) :: r,rmax_vsv,rmax_vsh,depthmax_depth
integer :: maxindex(1)
+ real(kind=CUSTOM_REAL) :: ratio_bulk,ratio_betav,ratio_betah,ratio_eta
! ------------------------------------------------------------------------
! sets maximum update in this depth range
logical,parameter :: use_depth_maximum = .true.
@@ -830,6 +841,35 @@
if( norm_eta_old < 1.e-22 ) call exit_mpi(myrank,'norm old gradient eta is zero')
endif
+ ! Powell, 1977: checks orthogonality between old and new gradients
+ ! gets length of ( gamma_(n-1)^T * gamma_n )
+ norm_bulk = sum( kernel_bulk_old * kernel_bulk )
+ norm_betav = sum( kernel_betav_old * kernel_betav )
+ norm_betah = sum( kernel_betah_old * kernel_betah )
+ norm_eta = sum( kernel_eta_old * kernel_eta )
+
+ call mpi_reduce(norm_bulk,norm_bulk_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_betav,norm_betav_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_betah,norm_betah_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_eta,norm_eta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+
+ ! ratio: ( g_n * g_n-1 ) / ( g_n-1 * g_n-1)
+ ratio_bulk = norm_bulk_sum / norm_bulk_old
+ ratio_betav = norm_betav_sum / norm_betav_old
+ ratio_betah = norm_betah_sum / norm_betah_old
+ ratio_eta = norm_eta_sum / norm_eta_old
+
+ ! if ratio > 0.2 (empirical threshold value), then on should restart with a steepest descent
+ if( myrank == 0 ) then
+ print*,'Powell ratio:'
+ print*,' bulk : ',ratio_bulk
+ print*,' betav: ',ratio_betav
+ print*,' betah: ',ratio_betah
+ print*,' eta : ',ratio_eta
+ print*
+ endif
+
+
! difference kernel/gradients
! length ( ( gamma_n - gamma_(n-1))^T * lambda_n )
norm_bulk = sum( (kernel_bulk - kernel_bulk_old) * kernel_bulk )
@@ -867,11 +907,33 @@
alpha_betah = norm_betah / norm_betah_old
alpha_eta = norm_eta / norm_eta_old
+ ! only if contribution is positive it will be considered, otherwise
+ ! we set it to zero so that it becomes a steepest descent update
+ if( alpha_bulk < 0.0 ) then
+ alpha_bulk = 0.0
+ endif
+ if( alpha_betav < 0.0 ) then
+ alpha_betav = 0.0
+ endif
+ if( alpha_betah < 0.0 ) then
+ alpha_betah = 0.0
+ endif
+ if( alpha_eta < 0.0 ) then
+ alpha_eta = 0.0
+ endif
+
else
! calculates only a single steplength applied to all
alpha_all = (norm_bulk + norm_betav + norm_betah + norm_eta) &
/ (norm_bulk_old + norm_betav_old + norm_betah_old + norm_eta_old)
- ! sets each steplength to same single one
+
+ ! only if contribution is positive it will be considered, otherwise
+ ! we set it to zero so that it becomes a steepest descent update
+ if( alpha_all < 0.0 ) then
+ alpha_all = 0.0
+ endif
+
+ ! sets each steplength to same single one
alpha_bulk = alpha_all
alpha_betav = alpha_all
alpha_betah = alpha_all
@@ -883,14 +945,14 @@
print*,' betav: ',alpha_betav
print*,' betah: ',alpha_betah
print*,' eta : ',alpha_eta
- print*
+ print*
endif
! broadcast values from rank 0 to all others
call mpi_bcast(alpha_bulk,1,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
call mpi_bcast(alpha_betav,1,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
call mpi_bcast(alpha_betah,1,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
call mpi_bcast(alpha_eta,1,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
-
+
! initializes kernel maximum
depthmax(:) = 0._CUSTOM_REAL
@@ -926,7 +988,7 @@
iglob = ibool(i,j,k,ispec)
r = sqrt( x(iglob)*x(iglob) + y(iglob)*y(iglob) + z(iglob)*z(iglob) )
- ! stores maximum kernel betav/betah value in this depth slice,
+ ! stores maximum kernel betav/betah value in this depth slice,
! since betav/betah are most likely dominating
if( r < R_top .and. r > R_bottom ) then
! kernel betav value
@@ -940,7 +1002,7 @@
if( depthmax(2) < max_vsh ) then
depthmax(2) = max_vsh
depthmax_radius(2) = r
- endif
+ endif
endif
endif
@@ -979,7 +1041,7 @@
iglob = ibool(i,j,k,ispec)
r = sqrt( x(iglob)*x(iglob) + y(iglob)*y(iglob) + z(iglob)*z(iglob) )
- ! stores maximum kernel betav/betah value in this depth slice,
+ ! stores maximum kernel betav/betah value in this depth slice,
! since betav/betah are most likely dominating
if( r < R_top .and. r > R_bottom ) then
! kernel betav value
@@ -993,7 +1055,7 @@
if( depthmax(2) < max_vsh ) then
depthmax(2) = max_vsh
depthmax_radius(2) = r
- endif
+ endif
endif
endif
@@ -1048,7 +1110,7 @@
depthmax(2) = max_vsh
depthmax_radius(1) = rmax_vsv
depthmax_radius(2) = rmax_vsh
-
+
max = maxval(depthmax)
maxindex = maxloc(depthmax)
depthmax_depth = depthmax_radius(maxindex(1))
@@ -1367,11 +1429,11 @@
integral_betah,integral_eta
real(kind=CUSTOM_REAL) :: volume_glob,volume_glob_sum
! root-mean square values
- real(kind=CUSTOM_REAL) :: rms_vpv,rms_vph,rms_vsv,rms_vsh,rms_eta,rms_rho
+ real(kind=CUSTOM_REAL) :: rms_vpv,rms_vph,rms_vsv,rms_vsh,rms_eta,rms_rho
real(kind=CUSTOM_REAL) :: rms_vpv_sum,rms_vph_sum,rms_vsv_sum,rms_vsh_sum, &
rms_eta_sum,rms_rho_sum
real(kind=CUSTOM_REAL) :: dvpv,dvph,dvsv,dvsh,deta,drho
-
+
! Gauss-Lobatto-Legendre points of integration and weights
double precision, dimension(NGLLX) :: xigll, wxgll
double precision, dimension(NGLLY) :: yigll, wygll
@@ -1453,9 +1515,9 @@
rms_vpv = 0._CUSTOM_REAL
rms_vph = 0._CUSTOM_REAL
rms_vsv = 0._CUSTOM_REAL
- rms_vsh = 0._CUSTOM_REAL
- rms_eta = 0._CUSTOM_REAL
- rms_rho = 0._CUSTOM_REAL
+ rms_vsh = 0._CUSTOM_REAL
+ rms_eta = 0._CUSTOM_REAL
+ rms_rho = 0._CUSTOM_REAL
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -1502,11 +1564,11 @@
call exit_MPI(myrank,'error NaN')
endif
- ! root-mean square
+ ! root-mean square
! integrates relative perturbations ( dv / v using logarithm ) squared
dvpv = log( model_vpv_new(i,j,k,ispec) / model_vpv(i,j,k,ispec) ) ! alphav
rms_vpv = rms_vpv + volumel * dvpv*dvpv
-
+
dvph = log( model_vph_new(i,j,k,ispec) / model_vph(i,j,k,ispec) ) ! alphah
rms_vph = rms_vph + volumel * dvph*dvph
@@ -1518,7 +1580,7 @@
deta = log( model_eta_new(i,j,k,ispec) / model_eta(i,j,k,ispec) ) ! eta
rms_eta = rms_eta + volumel * deta*deta
-
+
drho = log( model_rho_new(i,j,k,ispec) / model_rho(i,j,k,ispec) ) ! rho
rms_rho = rms_rho + volumel * drho*drho
@@ -1565,7 +1627,7 @@
print*
endif
- ! root-mean square
+ ! root-mean square
call mpi_reduce(rms_vpv,rms_vpv_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(rms_vph,rms_vph_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(rms_vsv,rms_vsv_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso_iso.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso_iso.f90 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso_iso.f90 2012-03-22 01:34:38 UTC (rev 19844)
@@ -87,7 +87,7 @@
! volume
real(kind=CUSTOM_REAL), dimension(NGLOB) :: x, y, z
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
- integer, dimension(NSPEC) :: idoubling
+ logical, dimension(NSPEC) :: ispec_is_tiso
! gradient vector norm ( v^T * v )
real(kind=CUSTOM_REAL) :: norm_a,norm_beta,norm_rho
@@ -325,6 +325,8 @@
use model_update_iso
implicit none
+ integer, dimension(:), allocatable :: idummy
+ logical, dimension(:), allocatable :: ldummy
! vpv model
write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vpv.bin'
@@ -416,6 +418,10 @@
print*
endif
+ ! allocates temporary array
+ allocate(idummy(NSPEC),ldummy(NSPEC),stat=ier)
+ if( ier /= 0 ) stop 'error allocating ldummy'
+
! global addressing
write(m_file,'(a,i6.6,a)') 'topo/proc',myrank,'_reg1_solver_data_2.bin'
open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
@@ -427,9 +433,14 @@
read(11) y(1:nglob)
read(11) z(1:nglob)
read(11) ibool(:,:,:,1:nspec)
- read(11) idoubling(1:nspec)
+ read(11) idummy(1:nspec) ! idoubling
+ read(11) ldummy(1:nspec) ! is_on_a_slice_edge
+ read(11) ispec_is_tiso(1:nspec)
close(11)
+ deallocate(idummy,ldummy)
+
+
end subroutine read_model
!
@@ -859,11 +870,11 @@
real(kind=CUSTOM_REAL) :: volume_glob,volume_glob_sum
! root-mean square values
- real(kind=CUSTOM_REAL) :: rms_vpv,rms_vph,rms_vsv,rms_vsh,rms_eta,rms_rho
+ real(kind=CUSTOM_REAL) :: rms_vpv,rms_vph,rms_vsv,rms_vsh,rms_eta,rms_rho
real(kind=CUSTOM_REAL) :: rms_vpv_sum,rms_vph_sum,rms_vsv_sum,rms_vsh_sum, &
rms_eta_sum,rms_rho_sum
real(kind=CUSTOM_REAL) :: dvpv,dvph,dvsv,dvsh,deta,drho
-
+
! Gauss-Lobatto-Legendre points of integration and weights
double precision, dimension(NGLLX) :: xigll, wxgll
double precision, dimension(NGLLY) :: yigll, wygll
@@ -943,9 +954,9 @@
rms_vpv = 0._CUSTOM_REAL
rms_vph = 0._CUSTOM_REAL
rms_vsv = 0._CUSTOM_REAL
- rms_vsh = 0._CUSTOM_REAL
- rms_eta = 0._CUSTOM_REAL
- rms_rho = 0._CUSTOM_REAL
+ rms_vsh = 0._CUSTOM_REAL
+ rms_eta = 0._CUSTOM_REAL
+ rms_rho = 0._CUSTOM_REAL
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -979,11 +990,11 @@
norm_beta = norm_beta + kernel_b(i,j,k,ispec)**2
norm_rho = norm_rho + kernel_rho(i,j,k,ispec)**2
- ! root-mean square
+ ! root-mean square
! integrates relative perturbations ( dv / v using logarithm ) squared
dvpv = log( model_vpv_new(i,j,k,ispec) / model_vpv(i,j,k,ispec) ) ! alphav
rms_vpv = rms_vpv + volumel * dvpv*dvpv
-
+
dvph = log( model_vph_new(i,j,k,ispec) / model_vph(i,j,k,ispec) ) ! alphah
rms_vph = rms_vph + volumel * dvph*dvph
@@ -995,7 +1006,7 @@
deta = log( model_eta_new(i,j,k,ispec) / model_eta(i,j,k,ispec) ) ! eta
rms_eta = rms_eta + volumel * deta*deta
-
+
drho = log( model_rho_new(i,j,k,ispec) / model_rho(i,j,k,ispec) ) ! rho
rms_rho = rms_rho + volumel * drho*drho
@@ -1038,7 +1049,7 @@
print*
endif
- ! root-mean square
+ ! root-mean square
call mpi_reduce(rms_vpv,rms_vpv_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(rms_vph,rms_vph_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(rms_vsv,rms_vsv_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/xcombine_vol_data
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/xcombine_vol_data 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/xcombine_vol_data 2012-03-22 01:34:38 UTC (rev 19844)
@@ -1 +0,0 @@
-link ../topo_input/combine_vol_data/xcombine_vol_data
\ No newline at end of file
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/scripts/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/scripts/README 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/scripts/README 2012-03-22 01:34:38 UTC (rev 19844)
@@ -2,14 +2,16 @@
/ADJOINT_TOMO/iterate_adj/pangu/scripts/README
-These scripts are used in setting up the forward simulations and the kernel simulations. They are included as a guide only, and will need to be modified for a particular setup.
+These scripts are used in setting up the forward simulations and the kernel simulations.
+They are included as a guide only, and will need to be modified for a particular setup.
setup_forward_dir.pl
setup_kernel_dir.pl
copy_kernel_dir.pl
copy_SEM_dir.pl
-In both cases, they require the directory Master_Dir, which is used to link run directories for all the requested event IDs. The organization is as follows:
+In both cases, they require the directory Master_Dir, which is used to link run directories
+for all the requested event IDs. The organization is as follows:
BASIN_KERNEL/
setup_kernel_dir.pl
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/smooth_sem_fun
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/smooth_sem_fun 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/smooth_sem_fun 2012-03-22 01:34:38 UTC (rev 19844)
@@ -1 +0,0 @@
-link src/smooth_sem_fun
\ No newline at end of file
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/xcombine_vol_data
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/xcombine_vol_data 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/xcombine_vol_data 2012-03-22 01:34:38 UTC (rev 19844)
@@ -1 +0,0 @@
-link ../topo_input/combine_vol_data/xcombine_vol_data
\ No newline at end of file
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth_update/smooth_sem_fun
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth_update/smooth_sem_fun 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth_update/smooth_sem_fun 2012-03-22 01:34:38 UTC (rev 19844)
@@ -1 +0,0 @@
-link ../smooth/smooth_sem_fun
\ No newline at end of file
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth_update/xcombine_vol_data
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth_update/xcombine_vol_data 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth_update/xcombine_vol_data 2012-03-22 01:34:38 UTC (rev 19844)
@@ -1 +0,0 @@
-link ../topo_input/combine_vol_data/xcombine_vol_data
\ No newline at end of file
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/kernels_run
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/kernels_run 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/kernels_run 2012-03-22 01:34:38 UTC (rev 19844)
@@ -1 +0,0 @@
-link ../KERNELS_MODELS/kernels_run
\ No newline at end of file
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/sum_kernels
===================================================================
(Binary files differ)
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/sum_kernels_globe.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/sum_kernels_globe.f90 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/sum_kernels_globe.f90 2012-03-22 01:34:38 UTC (rev 19844)
@@ -137,39 +137,53 @@
implicit none
+ include 'mpif.h'
include 'constants_globe.h'
+ include 'precision_globe.h'
include 'values_from_mesher_globe.h'
character(len=150) :: kernel_name,kernel_list(1000)
- integer :: nker,myrank
+ integer :: nker,myrank,ier
! local parameters
character(len=150) :: k_file
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
kernel_crust_mantle,total_kernel
+ real(kind=CUSTOM_REAL) :: norm,norm_sum
integer :: iker,ios
total_kernel = 0._CUSTOM_REAL
do iker = 1, nker
- if(myrank==0) then
- write(*,*) 'reading in event kernel for: ',trim(kernel_name)
- write(*,*) ' ',iker, ' out of ', nker
- endif
-
- write(k_file,'(a,i6.6,a)') 'INPUT_KERNELS/'//trim(kernel_list(iker)) &
- //'/proc',myrank,'_'//trim(kernel_name)//'.bin'
+ ! user output
+ if(myrank==0) then
+ write(*,*) 'reading in event kernel for: ',trim(kernel_name)
+ write(*,*) ' ',iker, ' out of ', nker
+ endif
- open(12,file=trim(k_file),status='old',form='unformatted',action='read',iostat=ios)
- if( ios /= 0 ) then
- write(*,*) ' kernel not found:',trim(k_file)
- cycle
- endif
+ ! sensitivity kernel
+ kernel_crust_mantle = 0._CUSTOM_REAL
+ write(k_file,'(a,i6.6,a)') 'INPUT_KERNELS/'//trim(kernel_list(iker)) &
+ //'/proc',myrank,'_'//trim(kernel_name)//'.bin'
- read(12) kernel_crust_mantle
- close(12)
+ open(12,file=trim(k_file),status='old',form='unformatted',action='read',iostat=ios)
+ if( ios /= 0 ) then
+ write(*,*) ' kernel not found:',trim(k_file)
+ cycle
+ endif
+ read(12) kernel_crust_mantle
+ close(12)
- total_kernel(:,:,:,1:NSPEC_CRUST_MANTLE_ADJOINT) = &
- total_kernel(:,:,:,1:NSPEC_CRUST_MANTLE_ADJOINT) &
+ ! outputs norm of kernel
+ norm = sum( kernel_crust_mantle * kernel_crust_mantle )
+ call mpi_reduce(norm,norm_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ print*,' norm kernel: ',sqrt(norm_sum)
+ print*
+ endif
+
+ ! sums all kernels from each event
+ total_kernel(:,:,:,1:NSPEC_CRUST_MANTLE_ADJOINT) = &
+ total_kernel(:,:,:,1:NSPEC_CRUST_MANTLE_ADJOINT) &
+ kernel_crust_mantle(:,:,:,1:NSPEC_CRUST_MANTLE_ADJOINT)
enddo
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/sum_preconditioned_kernels_globe.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/sum_preconditioned_kernels_globe.f90 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/sum_preconditioned_kernels_globe.f90 2012-03-22 01:34:38 UTC (rev 19844)
@@ -160,9 +160,9 @@
logical, parameter :: USE_SOURCE_MASK = .false.
!----------------------------------------------------------------------------------------
-
+ real(kind=CUSTOM_REAL) :: norm,norm_sum
character(len=150) :: kernel_name,kernel_list(1000)
- integer :: nker,myrank
+ integer :: nker,myrank,ier
! local parameters
character(len=150) :: k_file
@@ -189,6 +189,7 @@
! loops over all event kernels
total_kernel = 0._CUSTOM_REAL
do iker = 1, nker
+ ! user output
if(myrank==0) then
write(*,*) 'reading in event kernel for: ',trim(kernel_name)
write(*,*) ' and preconditioner: ','reg1_hess_kernel'
@@ -196,6 +197,7 @@
endif
! sensitivity kernel / frechet derivative
+ kernel_crust_mantle = 0._CUSTOM_REAL
write(k_file,'(a,i6.6,a)') 'INPUT_KERNELS/'//trim(kernel_list(iker)) &
//'/proc',myrank,'_'//trim(kernel_name)//'.bin'
@@ -207,7 +209,16 @@
read(12) kernel_crust_mantle
close(12)
+ ! outputs norm of kernel
+ norm = sum( kernel_crust_mantle * kernel_crust_mantle )
+ call mpi_reduce(norm,norm_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ print*,' norm kernel: ',sqrt(norm_sum)
+ endif
+
+
! approximate Hessian
+ hess_crust_mantle = 0._CUSTOM_REAL
write(k_file,'(a,i6.6,a)') 'INPUT_KERNELS/'//trim(kernel_list(iker)) &
//'/proc',myrank,'_reg1_hess_kernel.bin'
@@ -219,6 +230,14 @@
read(12) hess_crust_mantle
close(12)
+ ! outputs norm of preconditioner
+ norm = sum( hess_crust_mantle * hess_crust_mantle )
+ call mpi_reduce(norm,norm_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ print*,' norm preconditioner: ',sqrt(norm_sum)
+ print*
+ endif
+
! note: we take absolute values for hessian (as proposed by Yang)
hess_crust_mantle = abs(hess_crust_mantle)
@@ -325,6 +344,13 @@
call exit_mpi(myrank,'error hessian too small')
endif
+ ! user output
+ if( myrank == 0 ) then
+ print*
+ print*,'hessian maximum: ',maxh_all
+ print*
+ endif
+
! normalizes hessian
! since hessian has absolute values, this scales between [0,1]
hess_matrix = hess_matrix / maxh_all
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/xcombine_vol_data
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/xcombine_vol_data 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/xcombine_vol_data 2012-03-22 01:34:38 UTC (rev 19844)
@@ -1 +0,0 @@
-link ../topo_input/combine_vol_data/xcombine_vol_data
\ No newline at end of file
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/README 2012-03-22 01:03:41 UTC (rev 19843)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/README 2012-03-22 01:34:38 UTC (rev 19844)
@@ -8,14 +8,22 @@
compute_misfit.m
-After having run the measurement code, you can look at a set of misfit measurements using compute_misfit.m. In order to tabulate various summaries of the dataset (picks by station, event, network, etc), the program reads in a concatenated list of CMTSOLUTION files, and also a STATIONS list.
+After having run the measurement code, you can look at a set of misfit
+measurements using compute_misfit.m. In order to tabulate various summaries
+of the dataset (picks by station, event, network, etc), the program reads in
+a concatenated list of CMTSOLUTION files, and also a STATIONS list.
------------------
subspace_specfem.m
-After having run the measurement code AND the /subspace_hessian code (on the cluster), you can run subspace_specfem.m. The output files will be placed in the subdirectories of OUTPUT, which can be generated using make_dirs.pl. (The user should keep a copy of these files somewhere not on SVN.)
+After having run the measurement code AND the /subspace_hessian code (on the cluster),
+you can run subspace_specfem.m. The output files will be placed in the subdirectories
+of OUTPUT, which can be generated using make_dirs.pl.
+(The user should keep a copy of these files somewhere not on SVN.)
-Typically, one will generate, say, 100 model updates for 100 different regularization parameters, save the mu-vectors describing the linear combination of kernels describing the model update, and then generate the model updates on the cluster, using /subspace_update.
+Typically, one will generate, say, 100 model updates for 100 different regularization parameters,
+save the mu-vectors describing the linear combination of kernels describing the model update,
+and then generate the model updates on the cluster, using /subspace_update.
------------------
More information about the CIG-COMMITS
mailing list