[cig-commits] r19846 - in seismo/3D/ADJOINT_TOMO/iterate_adj/cluster: model_vp_vs model_vp_vs/src smooth/src sum_kernel/src
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Thu Mar 22 07:24:31 PDT 2012
Author: danielpeter
Date: 2012-03-22 07:24:30 -0700 (Thu, 22 Mar 2012)
New Revision: 19846
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/compile.bash
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/config.h
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/constants.h
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/precision.h
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/smooth_sem_fun.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/values_from_mesher.h
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/readme_globe
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/Makefile
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/Makefile_globe
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_iso.f90
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/smooth/src/Makefile
seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/Makefile
Log:
adds missing source files to smooth/src; updates global routines in model_vp_vs/src; separates compilation of global and local routines to Makefile and Makefile_globe
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/readme_globe
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/readme_globe 2012-03-22 04:17:38 UTC (rev 19845)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/readme_globe 2012-03-22 14:24:30 UTC (rev 19846)
@@ -40,8 +40,12 @@
takes TRANSVERSE ISOTROPIC model and TRANSVERSE ISOTROPIC kernels and
makes a steepest descent model update
+add_model_globe_tiso_cg:
+ takes TRANSVERSE ISOTROPIC model and TRANSVERSE ISOTROPIC kernels and
+ makes a conjugate-gradient model update
+
steps:
------
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/Makefile
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/Makefile 2012-03-22 04:17:38 UTC (rev 19845)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/Makefile 2012-03-22 14:24:30 UTC (rev 19846)
@@ -1,32 +1,12 @@
FFLAGS=-O3 -assume byterecl
-all: add_model_globe_iso \
- add_model_globe_tiso \
- add_model_globe_tiso_cg \
- add_model_globe_tiso_srvm \
- add_model_globe_tiso_iso
+all: add_model
-add_model_globe_iso: add_model_globe_iso.f90
- mpif90 -o add_model_globe_iso $(FFLAGS) add_model_globe_iso.f90 gll_library.f90 exit_mpi.f90
+add_model: add_model.f90
+ mpif90 -o add_model $(FFLAGS) add_model.f90 gll_library.f90 exit_mpi.f90
-add_model_globe_tiso: add_model_globe_tiso.f90
- mpif90 -o add_model_globe_tiso $(FFLAGS) add_model_globe_tiso.f90 gll_library.f90 exit_mpi.f90
-add_model_globe_tiso_iso: add_model_globe_tiso_iso.f90
- mpif90 -o add_model_globe_tiso_iso $(FFLAGS) add_model_globe_tiso_iso.f90 gll_library.f90 exit_mpi.f90
-
-add_model_globe_tiso_cg: add_model_globe_tiso_cg.f90
- mpif90 -o add_model_globe_tiso_cg $(FFLAGS) add_model_globe_tiso_cg.f90 gll_library.f90 exit_mpi.f90
-
-add_model_globe_tiso_srvm: add_model_globe_tiso_srvm.f90
- mpif90 -o add_model_globe_tiso_srvm $(FFLAGS) add_model_globe_tiso_srvm.f90 gll_library.f90 exit_mpi.f90
-
-
clean:
- rm -f add_model add_model_globe_iso \
- add_model_globe_tiso \
- add_model_globe_tiso_cg \
- add_model_globe_tiso_srvm \
- *.o *.mod
+ rm -f add_model *.o *.mod
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/Makefile_globe
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/Makefile_globe 2012-03-22 04:17:38 UTC (rev 19845)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/Makefile_globe 2012-03-22 14:24:30 UTC (rev 19846)
@@ -1,23 +1,26 @@
-FC=mpif90
FFLAGS=-O3 -assume byterecl
all: add_model_globe_iso \
add_model_globe_tiso \
+ add_model_globe_tiso_cg \
add_model_globe_tiso_iso
add_model_globe_iso: add_model_globe_iso.f90
- $(FC) -o add_model_globe_iso $(FFLAGS) add_model_globe_iso.f90 gll_library.f90 exit_mpi.f90
+ mpif90 -o add_model_globe_iso $(FFLAGS) add_model_globe_iso.f90 gll_library.f90 exit_mpi.f90
add_model_globe_tiso: add_model_globe_tiso.f90
- $(FC) -o add_model_globe_tiso $(FFLAGS) add_model_globe_tiso.f90 gll_library.f90 exit_mpi.f90
+ mpif90 -o add_model_globe_tiso $(FFLAGS) add_model_globe_tiso.f90 gll_library.f90 exit_mpi.f90
add_model_globe_tiso_iso: add_model_globe_tiso_iso.f90
- $(FC) -o add_model_globe_tiso_iso $(FFLAGS) add_model_globe_tiso_iso.f90 gll_library.f90 exit_mpi.f90
+ mpif90 -o add_model_globe_tiso_iso $(FFLAGS) add_model_globe_tiso_iso.f90 gll_library.f90 exit_mpi.f90
+add_model_globe_tiso_cg: add_model_globe_tiso_cg.f90
+ mpif90 -o add_model_globe_tiso_cg $(FFLAGS) add_model_globe_tiso_cg.f90 gll_library.f90 exit_mpi.f90
clean:
rm -f add_model_globe_iso \
- add_model_globe_tiso \
- add_model_globe_tiso_iso \
- *.o *.mod
\ No newline at end of file
+ add_model_globe_tiso \
+ add_model_globe_tiso_cg \
+ add_model_globe_tiso_iso \
+ *.o *.mod
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_iso.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_iso.f90 2012-03-22 04:17:38 UTC (rev 19845)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_iso.f90 2012-03-22 14:24:30 UTC (rev 19846)
@@ -37,10 +37,8 @@
!
! USAGE: e.g. ./add_model_globe_iso 0.3
-program add_model
+module model_update_iso
- implicit none
-
include 'mpif.h'
include 'constants_globe.h'
include 'precision_globe.h'
@@ -66,53 +64,149 @@
integer, parameter :: NSPEC = NSPEC_CRUST_MANTLE
integer, parameter :: NGLOB = NGLOB_CRUST_MANTLE
- character(len=150) :: sline, m_file, fname
+ ! isotropic model files
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ model_vp,model_vs,model_rho
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ model_vp_new,model_vs_new,model_rho_new
- ! model values, model updates and kernels
+ ! model updates
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
- model_vp, model_vs, model_rho, &
- model_dA, model_dB, model_dR, &
- kernel_a, kernel_b, kernel_rho, total_model
+ model_dA,model_dB,model_dR
+ ! kernels
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ kernel_a,kernel_b,kernel_rho
+
+ ! volume
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: x, y, z
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
+
+ ! gradient vector norm ( v^T * v )
+ real(kind=CUSTOM_REAL) :: norm_bulk,norm_beta,norm_rho
+ real(kind=CUSTOM_REAL) :: norm_bulk_sum,norm_beta_sum,norm_rho_sum
+
! steepest descent lengths
real(kind=CUSTOM_REAL) :: step_fac,step_length
+ ! statistics
+ real(kind=CUSTOM_REAL) :: min_vp,min_vs,max_vp,max_vs, &
+ min_rho,max_rho,minmax(4),vs_sum,vp_sum,rho_sum
+
+ real(kind=CUSTOM_REAL) :: beta1,beta0,rho1,rho0,alpha1,alpha0
+ real(kind=CUSTOM_REAL) :: dbetaiso,dA
+
integer :: nfile, myrank, sizeprocs, ier
integer :: i, j, k,ispec, iglob, ishell, n, it, j1, ib, npts_sem, ios
+ character(len=150) :: sline, m_file, fname
- ! statistics
- real(kind=CUSTOM_REAL) :: min_vp,min_vs,max_vp,max_vs, &
- min_rho,max_rho,max,minmax(4),vs_sum,vp_sum,rho_sum
- character(len=150) :: s_step_fac
+end module model_update_iso
- ! volume
- real(kind=CUSTOM_REAL), dimension(NGLOB) :: x, y, z
- integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: volume_spec
- real(kind=CUSTOM_REAL), dimension(NGLOB) :: volume_glob,volume_glob_sum
+!
+!-------------------------------------------------------------------------------------------------
+!
- ! Gauss-Lobatto-Legendre points of integration and weights
- double precision, dimension(NGLLX) :: xigll, wxgll
- double precision, dimension(NGLLY) :: yigll, wygll
- double precision, dimension(NGLLZ) :: zigll, wzgll
+program add_model
- ! array with all the weights in the cube
- double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+ use model_update_iso
- ! jacobian
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: jacobian
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
- real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl, &
- jacobianl,volumel
+ implicit none
- ! calculates volume associated with (global) point
- logical,parameter:: VOLUME_GLOBALPOINT = .true.
- real(kind=CUSTOM_REAL) :: kernel_integral_alpha,kernel_integral_beta,kernel_integral_rho
- real(kind=CUSTOM_REAL) :: integral_alpha_sum,integral_beta_sum,integral_rho_sum
-
! ============ program starts here =====================
+ ! initializes arrays
+ call initialize()
+
+ ! reads in parameters needed
+ call read_parameters()
+
+ ! reads in current transverse isotropic model files: vpv.. & vsv.. & eta & rho
+ call read_model()
+
+ ! reads in smoothed kernels: bulk, beta, rho
+ call read_kernels()
+
+ ! calculates gradient
+ ! steepest descent method
+ call get_gradient()
+
+
+ ! computes new model values for alpha, beta and rho
+ ! and stores new model files
+
+ ! model update:
+ ! isotropic update everywhere
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ ! initial model values
+ beta0 = model_vs(i,j,k,ispec)
+ rho0 = model_rho(i,j,k,ispec)
+ alpha0 = model_vp(i,j,k,ispec)
+
+ beta1 = 0._CUSTOM_REAL
+ rho1 = 0._CUSTOM_REAL
+ alpha1 = 0._CUSTOM_REAL
+
+ ! isotropic model update
+
+ ! shear values
+ dbetaiso = model_dB(i,j,k,ispec)
+ beta1 = beta0 * exp( dbetaiso )
+
+ ! density
+ rho1 = rho0 * exp( model_dR(i,j,k,ispec) )
+
+ ! alpha values
+ dA = model_dA(i,j,k,ispec)
+ if( USE_ALPHA_BETA_RHO ) then
+ ! new vp values use alpha model update
+ alpha1 = alpha0 * exp( dA )
+ else
+ ! new vp values use bulk model update:
+ ! this is based on vp_new = sqrt( bulk_new**2 + 4/3 vs_new**2 )
+ alpha1 = sqrt( alpha0**2 * exp(2.0*dA) + FOUR_THIRDS * beta0**2 * ( &
+ exp(2.0*dbetaiso) - exp(2.0*dA) ) )
+ endif
+
+ ! stores new model values
+ model_vp_new(i,j,k,ispec) = alpha1
+ model_vs_new(i,j,k,ispec) = beta1
+ model_rho_new(i,j,k,ispec) = rho1
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! stores new model in files
+ call store_new_model()
+
+ ! stores relative model perturbations
+ call store_perturbations()
+
+ ! computes volume element associated with points, calculates kernel integral for statistics
+ call compute_volume()
+
+ ! stop all the MPI processes, and exit
+ call MPI_FINALIZE(ier)
+
+end program add_model
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine initialize()
+
+! initializes arrays
+
+ use model_update_iso
+ implicit none
+
! initialize the MPI communicator and start the NPROCTOT MPI processes
call MPI_INIT(ier)
call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
@@ -126,13 +220,29 @@
model_vp = 0.0_CUSTOM_REAL
model_vs = 0.0_CUSTOM_REAL
model_rho = 0.0_CUSTOM_REAL
+
model_dA = 0.0_CUSTOM_REAL
model_dB = 0.0_CUSTOM_REAL
model_dR = 0.0_CUSTOM_REAL
+
kernel_a = 0.0_CUSTOM_REAL
kernel_b = 0.0_CUSTOM_REAL
kernel_rho = 0.0_CUSTOM_REAL
+end subroutine initialize
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine read_parameters()
+
+! reads in parameters needed
+
+ use model_update_iso
+ implicit none
+ character(len=150) :: s_step_fac
+
!--------------------------------------------------------
! subjective step length to multiply to the gradient
@@ -146,10 +256,11 @@
! read in parameter information
read(s_step_fac,*) step_fac
- if( abs(step_fac) < 1.e-10) then
- print*,'error: step factor ',step_fac
- call exit_MPI(myrank,'error step factor')
- endif
+
+ !if( abs(step_fac) < 1.e-10) then
+ ! print*,'error: step factor ',step_fac
+ ! call exit_MPI(myrank,'error step factor')
+ !endif
if (myrank == 0) then
print*,'defaults'
@@ -170,12 +281,25 @@
print*
endif
- open(20,file='step_fac',status='unknown',action='write')
- write(20,'(1e24.12)') step_fac
- close(20)
+ !open(20,file='step_fac',status='unknown',action='write')
+ !write(20,'(1e24.12)') step_fac
+ !close(20)
endif
- !-----------------------------------------------------
+end subroutine read_parameters
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine read_model()
+
+! reads in current transverse isotropic model: vpv.. & vsv.. & eta & rho
+
+ use model_update_iso
+ implicit none
+
+
! reads in current vp & vs & rho model
! vp model
write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vp.bin'
@@ -225,6 +349,34 @@
print*
endif
+ ! 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)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi(myrank,'file not found')
+ endif
+ read(11) x(1:nglob)
+ read(11) y(1:nglob)
+ read(11) z(1:nglob)
+ read(11) ibool(:,:,:,1:nspec)
+ close(11)
+
+end subroutine read_model
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine read_kernels()
+
+! reads in smoothed kernels: bulk, betav, betah, eta
+
+ use model_update_iso
+ implicit none
+
+
+
! reads in smoothed (& summed) event kernel
if( USE_ALPHA_BETA_RHO ) then
! reads in alpha kernel
@@ -296,134 +448,34 @@
print*
endif
+end subroutine read_kernels
- !--------------------------------------------------------
+!
+!-------------------------------------------------------------------------------------------------
+!
- if( VOLUME_GLOBALPOINT ) then
- ! computes volume element associated with points
- ! GLL points
- wgll_cube = 0.0
- call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
- call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
- call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- wgll_cube(i,j,k) = wxgll(i)*wygll(j)*wzgll(k)
- enddo
- enddo
- enddo
+subroutine get_gradient()
- ! 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)
- if( ier /= 0 ) then
- print*,'error opening: ',trim(m_file)
- call exit_mpi(myrank,'file not found')
- endif
- read(11) x(1:nglob)
- read(11) y(1:nglob)
- read(11) z(1:nglob)
- read(11) ibool(:,:,:,1:nspec)
- close(11)
+! calculates gradient by steepest descent method
- ! builds jacobian
- write(m_file,'(a,i6.6,a)') 'topo/proc',myrank,'_reg1_solver_data_1.bin'
- open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
- if( ier /= 0 ) then
- print*,'error opening: ',trim(m_file)
- call exit_mpi(myrank,'file not found')
- endif
- read(11) xix
- read(11) xiy
- read(11) xiz
- read(11) etax
- read(11) etay
- read(11) etaz
- read(11) gammax
- read(11) gammay
- read(11) gammaz
- close(11)
+ use model_update_iso
+ implicit none
+ ! local parameters
+ real(kind=CUSTOM_REAL):: r,max,depth_max
- jacobian = 0.0
- do ispec = 1, NSPEC
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- ! gets derivatives of ux, uy and uz with respect to x, y and z
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- ! computes the jacobian
- jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
- - xiyl*(etaxl*gammazl-etazl*gammaxl) &
- + xizl*(etaxl*gammayl-etayl*gammaxl))
- jacobian(i,j,k,ispec) = jacobianl
+ ! ------------------------------------------------------------------------
- enddo
- enddo
- enddo
- enddo
+ ! sets maximum update in this depth range
+ logical,parameter :: use_depth_maximum = .true.
+ ! normalized radii
+ real(kind=CUSTOM_REAL),parameter :: R_top = (6371.0 - 50.0 ) / R_EARTH_KM ! shallow depth
+ real(kind=CUSTOM_REAL),parameter :: R_bottom = (6371.0 - 100.0 ) / R_EARTH_KM ! deep depth
- ! volume associated with global point
- volume_glob = 0.0
- kernel_integral_alpha = 0._CUSTOM_REAL
- kernel_integral_beta = 0._CUSTOM_REAL
- kernel_integral_rho = 0._CUSTOM_REAL
- do ispec = 1, NSPEC
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- iglob = ibool(i,j,k,ispec)
- if( iglob == 0 ) then
- print*,'iglob zero',i,j,k,ispec
- print*
- print*,'ibool:',ispec
- print*,ibool(:,:,:,ispec)
- print*
- call exit_MPI(myrank,'error ibool')
- endif
+ ! ------------------------------------------------------------------------
- ! volume associated with GLL point
- volumel = jacobian(i,j,k,ispec)*wgll_cube(i,j,k)
+ ! initializes kernel maximum
+ max = 0._CUSTOM_REAL
- ! kernel integration: for each element
- kernel_integral_alpha = kernel_integral_alpha &
- + volumel * kernel_a(i,j,k,ispec)
-
- kernel_integral_beta = kernel_integral_beta &
- + volumel * kernel_b(i,j,k,ispec)
-
- kernel_integral_rho = kernel_integral_rho &
- + volumel * kernel_rho(i,j,k,ispec)
-
- enddo
- enddo
- enddo
- enddo
-
- ! statistics
- ! kernel integration: for whole volume
- call mpi_reduce(kernel_integral_alpha,integral_alpha_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(kernel_integral_beta,integral_beta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(kernel_integral_rho,integral_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
-
- if( myrank == 0 ) then
- print*,'integral kernels:'
- print*,' beta : ',integral_beta_sum
- print*,' a : ',integral_alpha_sum
- print*,' rho : ',integral_rho_sum
- print*
- endif
-
- endif
-
! calculates gradient
do ispec = 1, NSPEC
do k = 1, NGLLZ
@@ -441,6 +493,23 @@
! for rho
model_dR(i,j,k,ispec) = - kernel_rho(i,j,k,ispec)
+ ! determines maximum kernel beta value within given radius
+ if( use_depth_maximum ) then
+ ! get radius of point
+ iglob = ibool(i,j,k,ispec)
+ r = sqrt( x(iglob)*x(iglob) + y(iglob)*y(iglob) + z(iglob)*z(iglob) )
+
+ ! stores maximum kernel betav value in this depth slice, since betav is most likely dominating
+ if( r < R_top .and. r > R_bottom ) then
+ ! shear kernel value
+ max_vs = abs( kernel_b(i,j,k,ispec) )
+ if( max < max_vs ) then
+ max = max_vs
+ depth_max = r
+ endif
+ endif
+ endif
+
enddo
enddo
enddo
@@ -458,21 +527,42 @@
if( myrank == 0 ) then
print*,'initial gradients:'
+ print*,' a min/max: ',min_vp,max_vp
print*,' beta min/max : ',min_vs,max_vs
- print*,' a min/max: ',min_vp,max_vp
print*,' rho min/max: ',min_rho,max_rho
print*
endif
+ ! determines maximum kernel betav value within given radius
+ if( use_depth_maximum ) then
+ ! maximum of all processes stored in max_vsv
+ call mpi_reduce(max,max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ max = max_vs
+ depth_max = 6371.0 *( 1.0 - depth_max )
+ endif
+
! master determines step length based on maximum gradient value (either vp or vs)
if( myrank == 0 ) then
- minmax(1) = abs(min_vs)
- minmax(2) = abs(max_vs)
+
+ ! determines maximum kernel betav value within given radius
+ if( use_depth_maximum ) then
+ print*,' using depth maximum between 50km and 100km: ',max
+ print*,' approximate depth maximum: ',depth_max
+ print*
+ else
+ ! maximum gradient values
+ minmax(1) = abs(min_vs)
+ minmax(2) = abs(max_vs)
+ minmax(3) = abs(min_vp)
+ minmax(4) = abs(max_vp)
- minmax(3) = abs(min_vp)
- minmax(4) = abs(max_vp)
-
- max = maxval(minmax)
+ ! maximum value of all kernel maxima
+ max = maxval(minmax)
+ print*,' using maximum: ',max
+ print*
+ endif
+
+ ! chooses step length such that it becomes the desired, given step factor as inputted
step_length = step_fac/max
print*,' step length : ',step_length,max
@@ -485,24 +575,24 @@
max_vp = sum( model_dA * model_dA )
max_vs = sum( model_dB * model_dB )
max_rho = sum( model_dR * model_dR )
- max_vp = sqrt(max_vp)
- max_vs = sqrt(max_vs)
- max_rho = sqrt(max_rho)
call mpi_reduce(max_vp,vp_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(max_vs,vs_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(max_rho,rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ max_vp = sqrt(vp_sum)
+ max_vs = sqrt(vs_sum)
+ max_rho = sqrt(rho_sum)
+
if( myrank == 0 ) then
- print*,' initial beta length : ',vs_sum
- print*,' initial a length: ',vp_sum
- print*,' initial rho length: ',rho_sum
+ print*,'norm model updates:'
+ print*,' initial a length: ',max_vp
+ print*,' initial beta length : ',max_vs
+ print*,' initial rho length: ',max_rho
print*
endif
- ! model updates:
! multiply model updates by a subjective factor that will change the step
-
model_dA(:,:,:,:) = step_length * model_dA(:,:,:,:)
model_dB(:,:,:,:) = step_length * model_dB(:,:,:,:)
model_dR(:,:,:,:) = step_length * model_dR(:,:,:,:)
@@ -520,142 +610,330 @@
if( myrank == 0 ) then
print*,'scaled gradients:'
+ print*,' a min/max: ',min_vp,max_vp
print*,' beta min/max : ',min_vs,max_vs
- print*,' a min/max: ',min_vp,max_vp
print*,' rho min/max: ',min_rho,max_rho
print*
endif
- !-----------------------------------------------------
- ! new model:
- ! computes new model values for alpha, beta and rho
- ! and stores new model files
+end subroutine get_gradient
- ! S wavespeed model
- total_model = 0._CUSTOM_REAL
- total_model = model_vs * exp( model_dB )
+!
+!-------------------------------------------------------------------------------------------------
+!
- call mpi_reduce(maxval(total_model),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(minval(total_model),min_vs,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+subroutine store_new_model()
- fname = 'vs_new'
+! file output for new model
+
+ use model_update_iso
+ implicit none
+
+ ! vp model
+ call mpi_reduce(maxval(model_vp_new),max_vp,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_vp_new),min_vp,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'vp_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
- write(12) total_model(:,:,:,1:nspec)
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) model_vp_new
close(12)
- ! Rho density model
- total_model = 0._CUSTOM_REAL
- total_model = model_rho * exp( model_dR )
-
- call mpi_reduce(maxval(total_model),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(minval(total_model),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
-
- fname = 'rho_new'
+ ! vs model
+ call mpi_reduce(maxval(model_vs_new),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_vs_new),min_vs,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'vs_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
- write(12) total_model(:,:,:,1:nspec)
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) model_vs_new
close(12)
- ! P wavespeed model
- total_model = 0._CUSTOM_REAL
- if( USE_ALPHA_BETA_RHO ) then
- ! new vp values use alpha model update
- total_model = model_vp * exp( model_dA )
- else
- ! new vp values use bulk model update:
- ! this is based on vp_new = sqrt( bulk_new**2 + 4/3 vs_new**2 )
- total_model = sqrt( model_vp**2 * exp(2.0*model_dA) + FOUR_THIRDS * model_vs**2 *( &
- exp(2.0*model_dB) - exp(2.0*model_dA) ) )
- endif
- call mpi_reduce(maxval(total_model),max_vp,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(minval(total_model),min_vp,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
-
- fname = 'vp_new'
+ ! rho model
+ call mpi_reduce(maxval(model_rho_new),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_rho_new),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'rho_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
- write(12) total_model(:,:,:,1:nspec)
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) model_rho_new
close(12)
-
if( myrank == 0 ) then
print*,'new models:'
+ print*,' vp min/max: ',min_vp,max_vp
print*,' vs min/max: ',min_vs,max_vs
- print*,' vp min/max: ',min_vp,max_vp
print*,' rho min/max: ',min_rho,max_rho
print*
endif
- ! P & S & rho model update
- ! stores the linear approximations of the model updates
+end subroutine store_new_model
- ! stores relative Vp model perturbations
- where( model_vp /= 0.0 ) total_model = ( total_model - model_vp) / model_vp
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine store_perturbations()
+
+! file output for new model
+
+ use model_update_iso
+ implicit none
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: total_model
+
+ ! vp relative perturbations
+ ! logarithmic perturbation: log( v_new) - log( v_old) = log( v_new / v_old )
+ total_model = 0.0_CUSTOM_REAL
+ where( model_vp /= 0.0 ) total_model = log( model_vp_new / model_vp)
+ ! or
+ ! linear approximation: (v_new - v_old) / v_old
+ !where( model_vp /= 0.0 ) total_model = ( model_vp_new - model_vp) / model_vp
+
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvpvp.bin'
open(12,file=trim(m_file),form='unformatted',action='write')
- write(12) total_model(:,:,:,1:nspec)
+ write(12) total_model
close(12)
call mpi_reduce(maxval(total_model),max_vp,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(total_model),min_vp,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
- ! relative S model perturbations
- total_model = 0._CUSTOM_REAL
+ ! vs relative perturbations
+ total_model = 0.0_CUSTOM_REAL
+ where( model_vs /= 0.0 ) total_model = log( model_vs_new / model_vs)
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvsvs.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vs,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! rho relative model perturbations
+ total_model = 0.0_CUSTOM_REAL
+ where( model_rho /= 0.0 ) total_model = log( model_rho_new / model_rho)
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_drhorho.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'relative update:'
+ print*,' dvp/vp min/max: ',min_vp,max_vp
+ print*,' dvs/vs min/max: ',min_vs,max_vs
+ print*,' drho/rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+end subroutine store_perturbations
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine compute_volume()
+
+! computes volume element associated with points
+
+ use model_update_iso
+ implicit none
+ ! jacobian
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: jacobian
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl, &
+ jacobianl,volumel
+
+ ! integration values
+ real(kind=CUSTOM_REAL) :: kernel_integral_alpha,kernel_integral_beta,kernel_integral_rho
+ real(kind=CUSTOM_REAL) :: integral_alpha_sum,integral_beta_sum,integral_rho_sum
+
+ real(kind=CUSTOM_REAL) :: volume_glob,volume_glob_sum
+ ! root-mean square values
+ real(kind=CUSTOM_REAL) :: rms_vp,rms_vs,rms_rho
+ real(kind=CUSTOM_REAL) :: rms_vp_sum,rms_vs_sum,rms_rho_sum
+ real(kind=CUSTOM_REAL) :: dvp,dvs,drho
+
+ ! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: xigll, wxgll
+ double precision, dimension(NGLLY) :: yigll, wygll
+ double precision, dimension(NGLLZ) :: zigll, wzgll
+ ! array with all the weights in the cube
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+
+ ! GLL points
+ wgll_cube = 0.0
+ call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ wgll_cube(i,j,k) = wxgll(i)*wygll(j)*wzgll(k)
+ enddo
+ enddo
+ enddo
+
+ ! builds jacobian
+ write(m_file,'(a,i6.6,a)') 'topo/proc',myrank,'_reg1_solver_data_1.bin'
+ open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi(myrank,'file not found')
+ endif
+ read(11) xix
+ read(11) xiy
+ read(11) xiz
+ read(11) etax
+ read(11) etay
+ read(11) etaz
+ read(11) gammax
+ read(11) gammay
+ read(11) gammaz
+ close(11)
+
+ jacobian = 0.0
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
- if ( model_vs(i,j,k,ispec) > 1.e-10 ) then
- total_model(i,j,k,ispec) = ( model_vs(i,j,k,ispec)*exp(model_dB(i,j,k,ispec)) &
- - model_vs(i,j,k,ispec) ) / model_vs(i,j,k,ispec)
- endif
+ ! gets derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ ! computes the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+ jacobian(i,j,k,ispec) = jacobianl
+
enddo
enddo
enddo
enddo
- write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvsvs.bin'
- open(12,file=trim(m_file),form='unformatted',action='write')
- write(12) total_model(:,:,:,1:nspec)
- close(12)
- call mpi_reduce(maxval(total_model),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(minval(total_model),min_vs,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
- ! relative rho model perturbations
- total_model = 0._CUSTOM_REAL
+ ! volume associated with global point
+ volume_glob = 0._CUSTOM_REAL
+ kernel_integral_alpha = 0._CUSTOM_REAL
+ kernel_integral_beta = 0._CUSTOM_REAL
+ kernel_integral_rho = 0._CUSTOM_REAL
+ norm_bulk = 0._CUSTOM_REAL
+ norm_beta = 0._CUSTOM_REAL
+ norm_rho = 0._CUSTOM_REAL
+ rms_vp = 0._CUSTOM_REAL
+ rms_vs = 0._CUSTOM_REAL
+ rms_rho = 0._CUSTOM_REAL
+
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
- if ( model_rho(i,j,k,ispec) > 1.e-10 ) then
- total_model(i,j,k,ispec) = ( model_rho(i,j,k,ispec)*exp(model_dR(i,j,k,ispec)) &
- - model_rho(i,j,k,ispec) ) / model_rho(i,j,k,ispec)
+ iglob = ibool(i,j,k,ispec)
+ if( iglob == 0 ) then
+ print*,'iglob zero',i,j,k,ispec
+ print*
+ print*,'ibool:',ispec
+ print*,ibool(:,:,:,ispec)
+ print*
+ call exit_MPI(myrank,'error ibool')
endif
+
+ ! volume associated with GLL point
+ volumel = jacobian(i,j,k,ispec)*wgll_cube(i,j,k)
+ volume_glob = volume_glob + volumel
+
+ ! kernel integration: for each element
+ kernel_integral_alpha = kernel_integral_alpha &
+ + volumel * kernel_a(i,j,k,ispec)
+
+ kernel_integral_beta = kernel_integral_beta &
+ + volumel * kernel_b(i,j,k,ispec)
+
+ kernel_integral_rho = kernel_integral_rho &
+ + volumel * kernel_rho(i,j,k,ispec)
+
+ ! gradient vector norm sqrt( v^T * v )
+ norm_bulk = norm_bulk + kernel_a(i,j,k,ispec)**2
+ norm_beta = norm_beta + kernel_b(i,j,k,ispec)**2
+ norm_rho = norm_rho + kernel_rho(i,j,k,ispec)**2
+
+ ! checks number
+ if( isNaN(kernel_integral_alpha) ) then
+ print*,'error NaN: ',kernel_integral_alpha
+ print*,'rank:',myrank
+ print*,'i,j,k,ispec:',i,j,k,ispec
+ print*,'volumel: ',volumel,'kernel_bulk:',kernel_a(i,j,k,ispec)
+ call exit_MPI(myrank,'error NaN')
+ endif
+
+ ! root-mean square
+ ! integrates relative perturbations ( dv / v using logarithm ) squared
+ dvp = log( model_vp_new(i,j,k,ispec) / model_vp(i,j,k,ispec) ) ! alphav
+ rms_vp = rms_vp + volumel * dvp*dvp
+
+ dvs = log( model_vs_new(i,j,k,ispec) / model_vs(i,j,k,ispec) ) ! betav
+ rms_vs = rms_vs + volumel * dvs*dvs
+
+ drho = log( model_rho_new(i,j,k,ispec) / model_rho(i,j,k,ispec) ) ! rho
+ rms_rho = rms_rho + volumel * drho*drho
+
enddo
enddo
enddo
enddo
- write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_drhorho.bin'
- open(12,file=trim(m_file),form='unformatted',action='write')
- write(12) total_model(:,:,:,1:nspec)
- close(12)
- call mpi_reduce(maxval(total_model),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(minval(total_model),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ ! statistics
+ ! kernel integration: for whole volume
+ call mpi_reduce(kernel_integral_alpha,integral_alpha_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(kernel_integral_beta,integral_beta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(kernel_integral_rho,integral_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(volume_glob,volume_glob_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
if( myrank == 0 ) then
- print*,'relative Vs & Vp update:'
- print*,' dvs/vs min/max: ',min_vs,max_vs
- print*,' dvp/vp min/max: ',min_vp,max_vp
- print*,' drho/rho min/max: ',min_rho,max_rho
+ print*,'integral kernels:'
+ print*,' a : ',integral_alpha_sum
+ print*,' beta : ',integral_beta_sum
+ print*,' rho : ',integral_rho_sum
print*
+ print*,' total volume:',volume_glob_sum
+ print*
endif
+ ! norms: for whole volume
+ call mpi_reduce(norm_bulk,norm_bulk_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_beta,norm_beta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_rho,norm_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
- !-----------------------------------------------------
+ norm_bulk = sqrt(norm_bulk_sum)
+ norm_beta = sqrt(norm_beta_sum)
+ norm_rho = sqrt(norm_rho_sum)
- ! stop all the MPI processes, and exit
- call MPI_FINALIZE(ier)
+ if( myrank == 0 ) then
+ print*,'norm kernels:'
+ print*,' a : ',norm_bulk
+ print*,' beta : ',norm_beta
+ print*,' rho : ',norm_rho
+ print*
+ endif
-end program add_model
+ ! root-mean square
+ call mpi_reduce(rms_vp,rms_vp_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(rms_vs,rms_vs_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(rms_rho,rms_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ rms_vp = sqrt( rms_vp_sum / volume_glob_sum )
+ rms_vs = sqrt( rms_vs_sum / volume_glob_sum )
+ rms_rho = sqrt( rms_rho_sum / volume_glob_sum )
+ if( myrank == 0 ) then
+ print*,'root-mean square of perturbations:'
+ print*,' vp : ',rms_vp
+ print*,' vs : ',rms_vs
+ print*,' rho : ',rms_rho
+ print*
+ endif
+
+end subroutine compute_volume
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 04:17:38 UTC (rev 19845)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso.f90 2012-03-22 14:24:30 UTC (rev 19846)
@@ -91,7 +91,7 @@
real(kind=CUSTOM_REAL) :: min_vpv,min_vph,min_vsv,min_vsh, &
max_vpv,max_vph,max_vsv,max_vsh,min_eta,max_eta,min_bulk,max_bulk, &
- min_rho,max_rho,max,minmax(4)
+ min_rho,max_rho,minmax(4)
real(kind=CUSTOM_REAL) :: betav1,betah1,betav0,betah0,rho1,rho0, &
betaiso1,betaiso0,eta1,eta0,alphav1,alphav0,alphah1,alphah0
@@ -554,13 +554,16 @@
use model_update_tiso
implicit none
! local parameters
+ real(kind=CUSTOM_REAL):: r,max,depth_max
+
! ------------------------------------------------------------------------
+
! sets maximum update in this depth range
logical,parameter :: use_depth_maximum = .true.
! normalized radii
real(kind=CUSTOM_REAL),parameter :: R_top = (6371.0 - 50.0 ) / R_EARTH_KM ! shallow depth
real(kind=CUSTOM_REAL),parameter :: R_bottom = (6371.0 - 100.0 ) / R_EARTH_KM ! deep depth
- real(kind=CUSTOM_REAL):: r,depth_max
+
! ------------------------------------------------------------------------
! initializes kernel maximum
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 04:17:38 UTC (rev 19845)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso_cg.f90 2012-03-22 14:24:30 UTC (rev 19846)
@@ -794,7 +794,9 @@
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.
! normalized radii
@@ -859,14 +861,19 @@
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 ratio > 0.2 (empirical threshold value), then one should restart with a steepest descent
if( myrank == 0 ) then
- print*,'Powell ratio:'
+ print*,'Powell ratio: (> 0.2 then restart with steepest descent)'
print*,' bulk : ',ratio_bulk
print*,' betav: ',ratio_betav
print*,' betah: ',ratio_betah
print*,' eta : ',ratio_eta
print*
+ if( ratio_bulk > 0.2 .and. ratio_betav > 0.2 .and. ratio_betah > 0.2 &
+ .and. ratio_eta > 0.2 ) then
+ print*,' please consider doing a steepest descent instead cg...'
+ print*
+ endif
endif
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 04:17:38 UTC (rev 19845)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/model_vp_vs/src/add_model_globe_tiso_iso.f90 2012-03-22 14:24:30 UTC (rev 19846)
@@ -43,7 +43,7 @@
!
! USAGE: ./add_model_globe_tiso_iso 0.3
-module model_update_iso
+module model_update_tiso_iso
include 'mpif.h'
include 'constants_globe.h'
@@ -78,11 +78,11 @@
! model updates
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
- model_dA,model_dB,model_dR
+ model_dbulk,model_dbeta,model_drho
! kernels
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
- kernel_a,kernel_b,kernel_rho
+ kernel_bulk,kernel_beta,kernel_rho
! volume
real(kind=CUSTOM_REAL), dimension(NGLOB) :: x, y, z
@@ -90,15 +90,15 @@
logical, dimension(NSPEC) :: ispec_is_tiso
! gradient vector norm ( v^T * v )
- real(kind=CUSTOM_REAL) :: norm_a,norm_beta,norm_rho
- real(kind=CUSTOM_REAL) :: norm_a_sum,norm_beta_sum,norm_rho_sum
+ real(kind=CUSTOM_REAL) :: norm_bulk,norm_beta,norm_rho
+ real(kind=CUSTOM_REAL) :: norm_bulk_sum,norm_beta_sum,norm_rho_sum
! model update length
real(kind=CUSTOM_REAL) :: step_fac,step_length
real(kind=CUSTOM_REAL) :: min_vpv,min_vph,min_vsv,min_vsh, &
- max_vpv,max_vph,max_vsv,max_vsh,min_b,max_b,min_a,max_a, &
- min_rho,max_rho,min_eta,max_eta,max,minmax(4)
+ max_vpv,max_vph,max_vsv,max_vsh,min_beta,max_beta,min_bulk,max_bulk, &
+ min_rho,max_rho,min_eta,max_eta,minmax(4)
real(kind=CUSTOM_REAL) :: betav1,betah1,betav0,betah0,rho1,rho0, &
betaiso1,betaiso0,eta1,eta0,alphav1,alphav0,alphah1,alphah0
@@ -109,7 +109,7 @@
character(len=150) :: sline, m_file, fname
-end module model_update_iso
+end module model_update_tiso_iso
!
!-------------------------------------------------------------------------------------------------
@@ -117,7 +117,7 @@
program add_model
- use model_update_iso
+ use model_update_tiso_iso
implicit none
@@ -143,8 +143,7 @@
! (see also Carl's Latex notes)
! model update:
- ! transverse isotropic update only in layer Moho to 220 (where SPECFEM3D_GLOBE considers TISO)
- ! everywhere else uses an isotropic update
+ ! isotropic update everywhere
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -171,21 +170,29 @@
eta1 = eta0
! shear values
- dbetaiso = model_dB(i,j,k,ispec)
+ dbetaiso = model_dbeta(i,j,k,ispec)
betav1 = betav0 * exp( dbetaiso )
betah1 = betah0 * exp( dbetaiso )
! note: betah is probably not really used in isotropic layers
! (see SPECFEM3D_GLOBE/get_model.f90)
! density
- rho1 = rho0 * exp( model_dR(i,j,k,ispec) )
+ rho1 = rho0 * exp( model_drho(i,j,k,ispec) )
! alpha values
- dbulk = model_dA(i,j,k,ispec)
- alphav1 = sqrt( alphav0**2 * exp(2.0*dbulk) + FOUR_THIRDS * betav0**2 * ( &
+ dbulk = model_dbulk(i,j,k,ispec)
+ if( USE_ALPHA_BETA_RHO ) then
+ ! new vp values use alpha model update
+ alphav1 = alphav0 * exp( dbulk )
+ alphah1 = alphah0 * exp( dbulk )
+ else
+ ! new vp values use bulk model update:
+ ! this is based on vp_new = sqrt( bulk_new**2 + 4/3 vs_new**2 )
+ alphav1 = sqrt( alphav0**2 * exp(2.0*dbulk) + FOUR_THIRDS * betav0**2 * ( &
exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
- alphah1 = sqrt( alphah0**2 * exp(2.0*dbulk) + FOUR_THIRDS * betah0**2 * ( &
+ alphah1 = sqrt( alphah0**2 * exp(2.0*dbulk) + FOUR_THIRDS * betah0**2 * ( &
exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
+ endif
! note: alphah probably not used in SPECFEM3D_GLOBE
! stores new model values
@@ -223,7 +230,7 @@
! initializes arrays
- use model_update_iso
+ use model_update_tiso_iso
implicit none
! initialize the MPI communicator and start the NPROCTOT MPI processes
@@ -252,13 +259,13 @@
model_rho_new = 0.0_CUSTOM_REAL
! model updates
- model_dA = 0.0_CUSTOM_REAL
- model_dB = 0.0_CUSTOM_REAL
- model_dR = 0.0_CUSTOM_REAL
+ model_dbulk = 0.0_CUSTOM_REAL
+ model_dbeta = 0.0_CUSTOM_REAL
+ model_drho = 0.0_CUSTOM_REAL
! gradients
- kernel_a = 0.0_CUSTOM_REAL
- kernel_b = 0.0_CUSTOM_REAL
+ kernel_bulk = 0.0_CUSTOM_REAL
+ kernel_beta = 0.0_CUSTOM_REAL
kernel_rho = 0.0_CUSTOM_REAL
end subroutine initialize
@@ -271,7 +278,7 @@
! reads in parameters needed
- use model_update_iso
+ use model_update_tiso_iso
implicit none
character(len=150) :: s_step_fac
@@ -286,6 +293,7 @@
! read in parameter information
read(s_step_fac,*) step_fac
+
!if( abs(step_fac) < 1.e-10) then
! print*,'error: step factor ',step_fac
! call exit_MPI(myrank,'error step factor')
@@ -323,7 +331,7 @@
! reads in current transverse isotropic model: vpv.. & vsv.. & eta & rho
- use model_update_iso
+ use model_update_tiso_iso
implicit none
integer, dimension(:), allocatable :: idummy
logical, dimension(:), allocatable :: ldummy
@@ -451,7 +459,7 @@
! reads in smoothed kernels: bulk, betav, betah, eta
- use model_update_iso
+ use model_update_tiso_iso
implicit none
! bulk kernel
@@ -468,7 +476,7 @@
print*,'error opening: ',trim(m_file)
call exit_mpi(myrank,'file not found')
endif
- read(12) kernel_a(:,:,:,1:nspec)
+ read(12) kernel_bulk(:,:,:,1:nspec)
close(12)
! shear kernel
@@ -485,14 +493,14 @@
print*,'error opening: ',trim(m_file)
call exit_mpi(myrank,'file not found')
endif
- read(12) kernel_b(:,:,:,1:nspec)
+ read(12) kernel_beta(:,:,:,1:nspec)
close(12)
! rho kernel
if( USE_RHO_SCALING ) then
! uses scaling relation with shear perturbations
- kernel_rho(:,:,:,:) = RHO_SCALING * kernel_b(:,:,:,:)
+ kernel_rho(:,:,:,:) = RHO_SCALING * kernel_beta(:,:,:,:)
else
@@ -509,19 +517,19 @@
! statistics
- call mpi_reduce(minval(kernel_a),min_a,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(maxval(kernel_a),max_a,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(kernel_bulk),min_bulk,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(kernel_bulk),max_bulk,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(minval(kernel_b),min_b,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(maxval(kernel_b),max_b,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(kernel_beta),min_beta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(kernel_beta),max_beta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(kernel_rho),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
call mpi_reduce(maxval(kernel_rho),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
if( myrank == 0 ) then
print*,'initial kernels:'
- print*,' a min/max : ',min_a,max_a
- print*,' beta min/max: ',min_b,max_b
+ print*,' a min/max : ',min_bulk,max_bulk
+ print*,' beta min/max: ',min_beta,max_beta
print*,' rho min/max : ',min_rho,max_rho
print*
endif
@@ -537,9 +545,24 @@
! calculates gradient by steepest descent method
- use model_update_iso
+ use model_update_tiso_iso
implicit none
+ ! local parameters
+ real(kind=CUSTOM_REAL):: r,max,depth_max
+ ! ------------------------------------------------------------------------
+
+ ! sets maximum update in this depth range
+ logical,parameter :: use_depth_maximum = .true.
+ ! normalized radii
+ real(kind=CUSTOM_REAL),parameter :: R_top = (6371.0 - 50.0 ) / R_EARTH_KM ! shallow depth
+ real(kind=CUSTOM_REAL),parameter :: R_bottom = (6371.0 - 100.0 ) / R_EARTH_KM ! deep depth
+
+ ! ------------------------------------------------------------------------
+
+ ! initializes kernel maximum
+ max = 0._CUSTOM_REAL
+
! gradient in negative direction for steepest descent
do ispec = 1, NSPEC
do k = 1, NGLLZ
@@ -547,14 +570,31 @@
do i = 1, NGLLX
! for bulk
- model_dA(i,j,k,ispec) = - kernel_a(i,j,k,ispec)
+ model_dbulk(i,j,k,ispec) = - kernel_bulk(i,j,k,ispec)
! for shear
- model_dB(i,j,k,ispec) = - kernel_b(i,j,k,ispec)
+ model_dbeta(i,j,k,ispec) = - kernel_beta(i,j,k,ispec)
! for rho
- model_dR(i,j,k,ispec) = - kernel_rho(i,j,k,ispec)
+ model_drho(i,j,k,ispec) = - kernel_rho(i,j,k,ispec)
+ ! determines maximum kernel beta value within given radius
+ if( use_depth_maximum ) then
+ ! get radius of point
+ iglob = ibool(i,j,k,ispec)
+ r = sqrt( x(iglob)*x(iglob) + y(iglob)*y(iglob) + z(iglob)*z(iglob) )
+
+ ! stores maximum kernel betav value in this depth slice, since betav is most likely dominating
+ if( r < R_top .and. r > R_bottom ) then
+ ! shear kernel value
+ max_beta = abs( kernel_beta(i,j,k,ispec) )
+ if( max < max_beta ) then
+ max = max_beta
+ depth_max = r
+ endif
+ endif
+ endif
+
enddo
enddo
enddo
@@ -564,33 +604,53 @@
call store_kernel_updates()
! statistics
- call mpi_reduce(minval(model_dA),min_a,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(maxval(model_dA),max_a,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_dbulk),min_bulk,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dbulk),max_bulk,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(minval(model_dB),min_b,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(maxval(model_dB),max_b,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_dbeta),min_beta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dbeta),max_beta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(minval(model_dR),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(maxval(model_dR),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_drho),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_drho),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
if( myrank == 0 ) then
print*,'initial gradients:'
- print*,' a min/max : ',min_a,max_a
- print*,' beta min/max: ',min_b,max_b
+ print*,' a min/max : ',min_bulk,max_bulk
+ print*,' beta min/max: ',min_beta,max_beta
print*,' rho min/max : ',min_rho,max_rho
print*
endif
+ ! determines maximum kernel betav value within given radius
+ if( use_depth_maximum ) then
+ ! maximum of all processes stored in max_vsv
+ call mpi_reduce(max,max_beta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ max = max_beta
+ depth_max = 6371.0 *( 1.0 - depth_max )
+ endif
+
! determines step length based on maximum gradient value (either shear or bulk)
if( myrank == 0 ) then
- ! maximum gradient values
- minmax(1) = abs(min_b)
- minmax(2) = abs(max_b)
- minmax(3) = abs(min_a)
- minmax(4) = abs(max_a)
+ ! determines maximum kernel betav value within given radius
+ if( use_depth_maximum ) then
+ print*,' using depth maximum between 50km and 100km: ',max
+ print*,' approximate depth maximum: ',depth_max
+ print*
+ else
+ ! maximum gradient values
+ minmax(1) = abs(min_beta)
+ minmax(2) = abs(max_beta)
+ minmax(3) = abs(min_bulk)
+ minmax(4) = abs(max_bulk)
+
+ ! maximum value of all kernel maxima
+ max = maxval(minmax)
+ print*,' using maximum: ',max
+ print*
+ endif
+
! chooses step length such that it becomes the desired, given step factor as inputted
- max = maxval(minmax)
step_length = step_fac/max
print*,' step length : ',step_length,max
@@ -600,46 +660,46 @@
! gradient length sqrt( v^T * v )
- norm_a = sum( model_dA * model_dA )
- norm_beta = sum( model_dB * model_dB )
- norm_rho = sum( model_dR * model_dR )
+ norm_bulk = sum( model_dbulk * model_dbulk )
+ norm_beta = sum( model_dbeta * model_dbeta )
+ norm_rho = sum( model_drho * model_drho )
- call mpi_reduce(norm_a,norm_a_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_bulk,norm_bulk_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(norm_beta,norm_beta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(norm_rho,norm_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
- norm_a = sqrt(norm_a_sum)
+ norm_bulk = sqrt(norm_bulk_sum)
norm_beta = sqrt(norm_beta_sum)
norm_rho = sqrt(norm_rho_sum)
if( myrank == 0 ) then
print*,'norm model updates:'
- print*,' a : ',norm_a
+ print*,' a : ',norm_bulk
print*,' beta: ',norm_beta
print*,' rho : ',norm_rho
print*
endif
! multiply model updates by a subjective factor that will change the step
- model_dA(:,:,:,:) = step_length * model_dA(:,:,:,:)
- model_dB(:,:,:,:) = step_length * model_dB(:,:,:,:)
- model_dR(:,:,:,:) = step_length * model_dR(:,:,:,:)
+ model_dbulk(:,:,:,:) = step_length * model_dbulk(:,:,:,:)
+ model_dbeta(:,:,:,:) = step_length * model_dbeta(:,:,:,:)
+ model_drho(:,:,:,:) = step_length * model_drho(:,:,:,:)
! statistics
- call mpi_reduce(minval(model_dA),min_a,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(maxval(model_dA),max_a,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_dbulk),min_bulk,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dbulk),max_bulk,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(minval(model_dB),min_b,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(maxval(model_dB),max_b,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_dbeta),min_beta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dbeta),max_beta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(minval(model_dR),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(maxval(model_dR),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_drho),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_drho),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
if( myrank == 0 ) then
print*,'scaled gradients:'
- print*,' a min/max : ',min_a,max_a
- print*,' beta min/max: ',min_b,max_b
+ print*,' a min/max : ',min_bulk,max_bulk
+ print*,' beta min/max: ',min_beta,max_beta
print*,' rho min/max : ',min_rho,max_rho
print*
endif
@@ -654,26 +714,26 @@
! file output for new model
- use model_update_iso
+ use model_update_tiso_iso
implicit none
! kernel updates
- fname = 'dA'
+ fname = 'dbulk'
write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
- write(12) model_dA
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) model_dbulk
close(12)
- fname = 'dB'
+ fname = 'dbeta'
write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
- write(12) model_dB
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) model_dbeta
close(12)
- fname = 'dR'
+ fname = 'drho'
write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
- write(12) model_dR
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) model_drho
close(12)
end subroutine store_kernel_updates
@@ -686,7 +746,7 @@
! file output for new model
- use model_update_iso
+ use model_update_tiso_iso
implicit none
! vpv model
@@ -694,7 +754,7 @@
call mpi_reduce(minval(model_vpv_new),min_vpv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'vpv_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_vpv_new
close(12)
@@ -703,7 +763,7 @@
call mpi_reduce(minval(model_vph_new),min_vph,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'vph_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_vph_new
close(12)
@@ -712,7 +772,7 @@
call mpi_reduce(minval(model_vsv_new),min_vsv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'vsv_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_vsv_new
close(12)
@@ -721,7 +781,7 @@
call mpi_reduce(minval(model_vsh_new),min_vsh,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'vsh_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_vsh_new
close(12)
@@ -730,7 +790,7 @@
call mpi_reduce(minval(model_eta_new),min_eta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'eta_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_eta_new
close(12)
@@ -739,7 +799,7 @@
call mpi_reduce(minval(model_rho_new),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'rho_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_rho_new
close(12)
@@ -766,7 +826,7 @@
! file output for new model
- use model_update_iso
+ use model_update_tiso_iso
implicit none
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: total_model
@@ -856,7 +916,7 @@
! computes volume element associated with points
- use model_update_iso
+ use model_update_tiso_iso
implicit none
! jacobian
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: jacobian
@@ -865,8 +925,8 @@
real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl, &
jacobianl,volumel
! integration values
- real(kind=CUSTOM_REAL) :: integral_a_sum,integral_b_sum,integral_rho_sum
- real(kind=CUSTOM_REAL) :: kernel_integral_a,kernel_integral_b,kernel_integral_rho
+ real(kind=CUSTOM_REAL) :: integral_bulk_sum,integral_beta_sum,integral_rho_sum
+ real(kind=CUSTOM_REAL) :: integral_bulk,integral_beta,integral_rho
real(kind=CUSTOM_REAL) :: volume_glob,volume_glob_sum
! root-mean square values
@@ -945,10 +1005,10 @@
! volume associated with global point
volume_glob = 0._CUSTOM_REAL
- kernel_integral_a = 0._CUSTOM_REAL
- kernel_integral_b = 0._CUSTOM_REAL
- kernel_integral_rho = 0._CUSTOM_REAL
- norm_a = 0._CUSTOM_REAL
+ integral_bulk = 0._CUSTOM_REAL
+ integral_beta = 0._CUSTOM_REAL
+ integral_rho = 0._CUSTOM_REAL
+ norm_bulk = 0._CUSTOM_REAL
norm_beta = 0._CUSTOM_REAL
norm_rho = 0._CUSTOM_REAL
rms_vpv = 0._CUSTOM_REAL
@@ -957,6 +1017,7 @@
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
@@ -976,20 +1037,29 @@
volume_glob = volume_glob + volumel
! kernel integration: for each element
- kernel_integral_a = kernel_integral_a &
- + volumel * kernel_a(i,j,k,ispec)
+ integral_bulk = integral_bulk &
+ + volumel * kernel_bulk(i,j,k,ispec)
- kernel_integral_b = kernel_integral_b &
- + volumel * kernel_b(i,j,k,ispec)
+ integral_beta = integral_beta &
+ + volumel * kernel_beta(i,j,k,ispec)
- kernel_integral_rho = kernel_integral_rho &
+ integral_rho = integral_rho &
+ volumel * kernel_rho(i,j,k,ispec)
! gradient vector norm sqrt( v^T * v )
- norm_a = norm_a + kernel_a(i,j,k,ispec)**2
- norm_beta = norm_beta + kernel_b(i,j,k,ispec)**2
+ norm_bulk = norm_bulk + kernel_bulk(i,j,k,ispec)**2
+ norm_beta = norm_beta + kernel_beta(i,j,k,ispec)**2
norm_rho = norm_rho + kernel_rho(i,j,k,ispec)**2
+ ! checks number
+ if( isNaN(integral_bulk) ) then
+ print*,'error NaN: ',integral_bulk
+ print*,'rank:',myrank
+ print*,'i,j,k,ispec:',i,j,k,ispec
+ print*,'volumel: ',volumel,'kernel_bulk:',kernel_bulk(i,j,k,ispec)
+ call exit_MPI(myrank,'error NaN')
+ endif
+
! 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
@@ -1017,15 +1087,15 @@
! statistics
! kernel integration: for whole volume
- call mpi_reduce(kernel_integral_a,integral_a_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(kernel_integral_b,integral_b_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(kernel_integral_rho,integral_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(integral_bulk,integral_bulk_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(integral_beta,integral_beta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(integral_rho,integral_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(volume_glob,volume_glob_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
if( myrank == 0 ) then
print*,'integral kernels:'
- print*,' a : ',integral_a_sum
- print*,' beta : ',integral_b_sum
+ print*,' a : ',integral_bulk_sum
+ print*,' beta : ',integral_beta_sum
print*,' rho : ',integral_rho_sum
print*
print*,' total volume:',volume_glob_sum
@@ -1033,17 +1103,17 @@
endif
! norms: for whole volume
- call mpi_reduce(norm_a,norm_a_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_bulk,norm_bulk_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(norm_beta,norm_beta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(norm_rho,norm_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
- norm_a = sqrt(norm_a_sum)
+ norm_bulk = sqrt(norm_bulk_sum)
norm_beta = sqrt(norm_beta_sum)
norm_rho = sqrt(norm_rho_sum)
if( myrank == 0 ) then
print*,'norm kernels:'
- print*,' a : ',norm_a
+ print*,' a : ',norm_bulk
print*,' beta : ',norm_beta
print*,' rho : ',norm_rho
print*
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/Makefile
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/Makefile 2012-03-22 04:17:38 UTC (rev 19845)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/Makefile 2012-03-22 14:24:30 UTC (rev 19846)
@@ -1,12 +1,12 @@
FC=mpif90
FFLAGS=-O3 -assume byterecl #-g -traceback
-all: smooth_sem_globe
+all: smooth_sem_fun
-smooth_sem_globe: smooth_sem_globe.f90
- $(FC) -o smooth_sem_globe $(FFLAGS) smooth_sem_globe.f90 smooth_sub.f90 exit_mpi.f90 gll_library.f90
+smooth_sem_fun: smooth_sem_fun.f90
+ $(FC) -o smooth_sem_fun $(FFLAGS) smooth_sem_fun.f90 smooth_sub.f90 exit_mpi.f90 gll_library.f90
clean:
- rm -f smooth_sem_globe *.o
+ rm -f smooth_sem_fun *.o
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/compile.bash
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/compile.bash (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/compile.bash 2012-03-22 14:24:30 UTC (rev 19846)
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+flags="-O3 -implicitnone -warn argument_checking -warn unused -warn declarations -check nobounds"
+
+mpif90 $flags -o smooth_sem_fun smooth_sem_fun.f90 smooth_sub.f90 gll_library.f90 exit_mpi.f90
+
+#ifort $flags -o smooth_sem_fun smooth_specfem_function.f90 smooth_sub.f90 gll_library.f90 exit_mpi.f90
+
+rm *.o
+
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/compile.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/config.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/config.h (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/config.h 2012-03-22 14:24:30 UTC (rev 19846)
@@ -0,0 +1,35 @@
+/* config.h. Generated from config.h.in by configure. */
+/* config.h.in. Generated from configure.ac by autoheader. */
+
+/* Define to dummy `main' function (if any) required to link to the Fortran
+ libraries. */
+/* #undef FC_DUMMY_MAIN */
+
+/* Define if F77 and FC dummy `main' functions are identical. */
+/* #undef FC_DUMMY_MAIN_EQ_F77 */
+
+/* Define to a macro mangling the given C identifier (in lower and upper
+ case), which must not contain underscores, for linking with Fortran. */
+#define FC_FUNC(name,NAME) name ## _
+
+/* As FC_FUNC, but for C identifiers containing underscores. */
+#define FC_FUNC_(name,NAME) name ## _
+
+/* Define to alternate name for `main' routine that is called from a `main' in
+ the Fortran libraries. */
+/* #undef FC_MAIN */
+
+/* Define to the address where bug reports for this package should be sent. */
+#define PACKAGE_BUGREPORT "jtromp AT caltech.edu"
+
+/* Define to the full name of this package. */
+#define PACKAGE_NAME "Specfem 3D Basin"
+
+/* Define to the full name and version of this package. */
+#define PACKAGE_STRING "Specfem 3D Basin 1.4.0"
+
+/* Define to the one symbol short name of this package. */
+#define PACKAGE_TARNAME "Specfem3DBasin"
+
+/* Define to the version of this package. */
+#define PACKAGE_VERSION "1.4.0"
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/constants.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/constants.h (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/constants.h 2012-03-22 14:24:30 UTC (rev 19846)
@@ -0,0 +1,373 @@
+!=====================================================================
+!
+! S p e c f e m 3 D B a s i n V e r s i o n 1 . 4
+! --------------------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology July 2005
+!
+! A signed non-commercial agreement is required to use this program.
+! Please check http://www.gps.caltech.edu/research/jtromp for details.
+! Free for non-commercial academic research ONLY.
+! This program is distributed WITHOUT ANY WARRANTY whatsoever.
+! Do not redistribute this program without written permission.
+!
+!=====================================================================
+
+! constants.h. Generated from constants.h.in by configure.
+
+!
+! solver in single or double precision depending on the machine (4 or 8 bytes)
+!
+! ALSO CHANGE FILE precision.h ACCORDINGLY
+!
+ integer, parameter :: SIZE_REAL = 4
+ integer, parameter :: SIZE_DOUBLE = 8
+
+! set to SIZE_REAL to run in single precision
+! set to SIZE_DOUBLE to run in double precision (increases memory size by 2)
+ integer, parameter :: CUSTOM_REAL = SIZE_REAL
+
+!----------- parameters that can be changed by the user -----------
+
+! set to .false. if running on a Beowulf-type machine with local disks
+! set to .true. if running on a shared-memory machine with common file system
+! if running on a Beowulf, also modify name of nodes in filter_machine_file.f90
+ logical, parameter :: LOCAL_PATH_IS_ALSO_GLOBAL = .false.
+
+! apply heuristic rule to modify doubling regions to balance angles
+ logical, parameter :: APPLY_HEURISTIC_RULE = .true.
+
+! number of GLL points in each direction of an element (degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLY = NGLLX
+ integer, parameter :: NGLLZ = NGLLX
+
+! input, output and main MPI I/O files
+ integer, parameter :: ISTANDARD_OUTPUT = 6
+ integer, parameter :: IIN = 40,IOUT = 41
+! uncomment this to write messages to a text file
+ integer, parameter :: IMAIN = 42
+! uncomment this to write messages to the screen
+! integer, parameter :: IMAIN = ISTANDARD_OUTPUT
+! I/O unit for source and receiver vtk file
+ integer, parameter :: IOVTK = 98
+
+! minimum thickness in meters to include the effect of the oceans
+! to avoid taking into account spurious oscillations in topography model
+ double precision, parameter :: MINIMUM_THICKNESS_3D_OCEANS = 10.d0
+
+! min and max density in the model
+ double precision, parameter :: DENSITY_MAX = 3000.d0
+ double precision, parameter :: DENSITY_MIN = 2000.d0
+
+! density of sea water
+ real(kind=CUSTOM_REAL), parameter :: RHO_OCEANS = 1020.0
+
+! extend basin model below threshold and above topography to make sure
+! there is no small gap between interpolated maps and sediments
+ logical, parameter :: EXTEND_VOXET_BELOW_BASEMENT = .true.
+ logical, parameter :: EXTEND_VOXET_ABOVE_TOPO = .true.
+ double precision, parameter :: DISTMAX_ASSUME_SEDIMENTS = 210.d0
+ integer, parameter :: NCELLS_EXTEND = 8
+
+! depth at which we start to honor the basement interface
+ double precision, parameter :: Z_THRESHOLD_HONOR_BASEMENT = -4700.d0
+
+! flag to print the details of source location
+ logical, parameter :: SHOW_DETAILS_LOCATE_SOURCE = .false.
+
+! maximum length of station and network name for receivers
+ integer, parameter :: MAX_LENGTH_STATION_NAME = 32
+ integer, parameter :: MAX_LENGTH_NETWORK_NAME = 8
+
+! number of sources to be gathered by MPI_Gather
+ integer, parameter :: NGATHER_SOURCES = 10000
+
+! source decay rate
+ double precision, parameter :: SOURCE_DECAY_RATE = 1.628d0
+
+! ---------------------------------------------------------------------------------------
+! LQY -- Following 3 variables stays here temporarily,
+! we need to move them to Par_file at a proper time
+! ---------------------------------------------------------------------------------------
+! save moho mesh and compute Moho boundary kernels
+ logical, parameter :: SAVE_MOHO_MESH = .true.
+
+! number of steps to save the state variables in the forward simulation,
+! to be used in the backward reconstruction in the presence of attenuation
+ integer, parameter :: NSTEP_Q_SAVE = 200
+
+! the scratch disk to save the state variables saved in the forward
+! simulation, this can be a global scratch disk in case you run out of
+! space on the local scratch disk
+ character(len=150), parameter :: LOCAL_PATH_Q = '/ibrixfs1/scratch/lqy'
+
+!------------------------------------------------------
+!----------- do not modify anything below -------------
+!------------------------------------------------------
+
+! on some processors (e.g. Pentiums) it is necessary to suppress underflows
+! by using a small initial field instead of zero
+ logical, parameter :: FIX_UNDERFLOW_PROBLEM = .true.
+
+! some useful constants
+ double precision, parameter :: PI = 3.141592653589793d0
+ double precision, parameter :: TWO_PI = 2.d0 * PI
+
+! 3-D simulation
+ integer, parameter :: NDIM = 3
+
+! dimension of the boundaries of the slices
+ integer, parameter :: NDIM2D = 2
+
+! number of nodes for 2D and 3D shape functions for hexahedra
+! we use 8-node mesh bricks, which are more stable than 27-node elements
+ integer, parameter :: NGNOD = 8, NGNOD2D = 4
+
+! a few useful constants
+ double precision, parameter :: ZERO = 0.d0,ONE = 1.d0,TWO = 2.d0,HALF = 0.5d0
+
+ real(kind=CUSTOM_REAL), parameter :: &
+ ONE_THIRD = 1._CUSTOM_REAL/3._CUSTOM_REAL, &
+ FOUR_THIRDS = 4._CUSTOM_REAL/3._CUSTOM_REAL
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.d+30,TINYVAL = 1.d-9
+
+! very large real value declared independently of the machine
+ real(kind=CUSTOM_REAL), parameter :: HUGEVAL_SNGL = 1.e+30_CUSTOM_REAL
+
+! very large integer value
+ integer, parameter :: HUGEINT = 100000000
+
+! define flag for elements
+ integer, parameter :: IFLAG_ONE_LAYER_TOPOGRAPHY = 1
+ integer, parameter :: IFLAG_BASEMENT_TOPO = 2
+ integer, parameter :: IFLAG_16km_BASEMENT = 3
+ integer, parameter :: IFLAG_MOHO_16km = 4
+ integer, parameter :: IFLAG_HALFSPACE_MOHO = 5
+
+! define flag for regions for attenuation
+ integer, parameter :: NUM_REGIONS_ATTENUATION = 13
+
+ integer, parameter :: IATTENUATION_SEDIMENTS_40 = 1
+ integer, parameter :: IATTENUATION_SEDIMENTS_50 = 2
+ integer, parameter :: IATTENUATION_SEDIMENTS_60 = 3
+ integer, parameter :: IATTENUATION_SEDIMENTS_70 = 4
+ integer, parameter :: IATTENUATION_SEDIMENTS_80 = 5
+ integer, parameter :: IATTENUATION_SEDIMENTS_90 = 6
+ integer, parameter :: IATTENUATION_SEDIMENTS_100 = 7
+ integer, parameter :: IATTENUATION_SEDIMENTS_110 = 8
+ integer, parameter :: IATTENUATION_SEDIMENTS_120 = 9
+ integer, parameter :: IATTENUATION_SEDIMENTS_130 = 10
+ integer, parameter :: IATTENUATION_SEDIMENTS_140 = 11
+ integer, parameter :: IATTENUATION_SEDIMENTS_150 = 12
+ integer, parameter :: IATTENUATION_BEDROCK = 13
+
+! Olsen's constant for Q_mu = constant * v_s attenuation rule
+ real, parameter :: OLSEN_ATTENUATION_RATIO = 0.05
+
+! number of standard linear solids in parallel for attenuation
+ integer, parameter :: N_SLS = 3
+
+! flag for the four edges of each slice and for the bottom edge
+ integer, parameter :: XI_MIN = 1
+ integer, parameter :: XI_MAX = 2
+ integer, parameter :: ETA_MIN = 3
+ integer, parameter :: ETA_MAX = 4
+ integer, parameter :: BOTTOM = 5
+
+! number of points per surface element
+ integer, parameter :: NGLLSQUARE = NGLLX * NGLLY
+
+! number of points per spectral element
+ integer, parameter :: NGLLCUBE = NGLLX * NGLLY * NGLLZ
+
+! for vectorization of loops
+ integer, parameter :: NGLLSQUARE_NDIM = NGLLSQUARE * NDIM
+ integer, parameter :: NGLLCUBE_NDIM = NGLLCUBE * NDIM
+
+! flag for projection from latitude/longitude to UTM, and back
+ integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
+
+! smallest real number on the Pentium and the SGI = 1.1754944E-38
+! largest real number on the Pentium and the SGI = 3.4028235E+38
+! small negligible initial value to avoid very slow underflow trapping
+! but not too small to avoid trapping on velocity and acceleration in Newmark
+ real(kind=CUSTOM_REAL), parameter :: VERYSMALLVAL = 1.E-24_CUSTOM_REAL
+
+! displacement threshold above which we consider the code became unstable
+ real(kind=CUSTOM_REAL), parameter :: STABILITY_THRESHOLD = 1.E+25_CUSTOM_REAL
+
+! geometrical tolerance for boundary detection
+ double precision, parameter :: SMALLVAL = 0.00001d0
+
+! do not use tags for MPI messages, use dummy tag instead
+ integer, parameter :: itag = 0,itag2 = 0
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.d0,GAUSSBETA = 0.d0
+
+! number of lines per source in CMTSOLUTION file
+ integer, parameter :: NLINES_PER_CMTSOLUTION_SOURCE = 13
+
+! number of iterations to solve the system for xi and eta
+ integer, parameter :: NUM_ITER = 4
+
+! size of topography and bathymetry file for Southern California
+ integer, parameter :: NX_TOPO_SOCAL = 1401,NY_TOPO_SOCAL = 1001
+ double precision, parameter :: ORIG_LAT_TOPO_SOCAL = 32.d0
+ double precision, parameter :: ORIG_LONG_TOPO_SOCAL = -121.d0
+ double precision, parameter :: DEGREES_PER_CELL_TOPO_SOCAL = 5.d0 / 1000.d0
+ character(len=100), parameter :: TOPO_FILE_SOCAL = 'DATA/la_topography/topo_bathy_final.dat'
+
+! size of topography and bathymetry file for Lacq (France) gas field
+ integer, parameter :: NX_TOPO_LACQ = 499,NY_TOPO_LACQ = 401
+ double precision, parameter :: ORIG_LAT_TOPO_LACQ = 100000.d0
+ double precision, parameter :: ORIG_LONG_TOPO_LACQ = 340400.d0
+ double precision, parameter :: DEGREES_PER_CELL_TOPO_LACQ = 100.d0
+ character(len=100), parameter :: TOPO_FILE_LACQ = 'DATA/lacq_thomas/mnt_Lacq_Lambert_final_dimitri.dat'
+
+! size of Lupei Zhu's Moho map file for Southern California
+ integer, parameter :: NX_MOHO = 71,NY_MOHO = 51
+ double precision, parameter :: ORIG_LAT_MOHO = 32.d0
+ double precision, parameter :: ORIG_LONG_MOHO = -121.d0
+ double precision, parameter :: DEGREES_PER_CELL_MOHO = 0.1d0
+
+! size of basement map file
+ integer, parameter :: NX_BASEMENT = 161,NY_BASEMENT = 144
+ double precision, parameter :: ORIG_X_BASEMENT = 316000.
+ double precision, parameter :: ORIG_Y_BASEMENT = 3655000.
+ double precision, parameter :: SPACING_X_BASEMENT = 1000.
+ double precision, parameter :: SPACING_Y_BASEMENT = 1000.
+
+!
+! new Gocad Voxets Peter July 29, 2002 - high-res and medium-res blocks
+!
+
+! size of the medium-resolution Gocad voxet
+ integer, parameter :: NX_GOCAD_MR = 194, NY_GOCAD_MR = 196, NZ_GOCAD_MR = 100
+
+ double precision, parameter :: ORIG_X_GOCAD_MR = 283000.
+ double precision, parameter :: ORIG_Y_GOCAD_MR = 3655000.
+ double precision, parameter :: ORIG_Z_GOCAD_MR = -15000.
+
+ double precision, parameter :: SPACING_X_GOCAD_MR = 1000.
+ double precision, parameter :: SPACING_Y_GOCAD_MR = 1000.
+ double precision, parameter :: SPACING_Z_GOCAD_MR = 200.
+
+! maximum size of model for tapering of transition between Hauksson and MR
+ double precision, parameter :: END_X_GOCAD_MR = ORIG_X_GOCAD_MR + SPACING_X_GOCAD_MR * (NX_GOCAD_MR - 1)
+ double precision, parameter :: END_Y_GOCAD_MR = ORIG_Y_GOCAD_MR + SPACING_Y_GOCAD_MR * (NY_GOCAD_MR - 1)
+
+! size of the high-resolution Gocad voxet
+ integer, parameter :: NX_GOCAD_HR = 185, NY_GOCAD_HR = 196, NZ_GOCAD_HR = 100
+
+ double precision, parameter :: ORIG_X_GOCAD_HR = 371052.25
+ double precision, parameter :: ORIG_Y_GOCAD_HR = 3725250.
+ double precision, parameter :: ORIG_Z_GOCAD_HR = -9500.
+
+ double precision, parameter :: SPACING_X_GOCAD_HR = 250.
+ double precision, parameter :: SPACING_Y_GOCAD_HR = 250.
+ double precision, parameter :: SPACING_Z_GOCAD_HR = 100.
+
+! maximum size of model for tapering of transition between HR and MR
+ double precision, parameter :: END_X_GOCAD_HR = ORIG_X_GOCAD_HR + SPACING_X_GOCAD_HR * (NX_GOCAD_HR - 1)
+ double precision, parameter :: END_Y_GOCAD_HR = ORIG_Y_GOCAD_HR + SPACING_Y_GOCAD_HR * (NY_GOCAD_HR - 1)
+
+! implement smooth transition between Hauksson, HR and MR Gocad blocks
+ logical, parameter :: TAPER_GOCAD_TRANSITIONS = .true.
+
+! Salton Sea Gocad voxet
+ integer, parameter :: GOCAD_ST_NU = 638, GOCAD_ST_NV = 219, GOCAD_ST_NW = 76
+ double precision, parameter :: GOCAD_ST_O_X = 720844.0, GOCAD_ST_O_Y = 3401799.250, &
+ GOCAD_ST_O_Z = -6354.334
+ double precision, parameter :: GOCAD_ST_U_X = -209197.89, GOCAD_ST_U_Y = 320741.71
+ double precision, parameter :: GOCAD_ST_V_X = 109670.74, GOCAD_ST_V_Y = 71530.72
+ double precision, parameter :: GOCAD_ST_W_Z = 7666.334
+ double precision, parameter :: GOCAD_ST_NO_DATA_VALUE = -99999
+
+!
+!--- larger Hauksson model for entire So-Cal, 15 km resolution
+!
+
+!! Hauksson (2000)
+!! number of non-constant layers
+! integer, parameter :: NLAYERS_HAUKSSON = 9
+!! depth of layers
+! double precision, parameter :: Z_HAUKSSON_LAYER_1 = -1000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_2 = -4000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_3 = -6000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_4 = -10000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_5 = -15000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_6 = -17000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_7 = -22000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_8 = -31000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_9 = -33000.d0
+!
+! integer, parameter :: NGRID_NEW_HAUKSSON = 201
+!
+!! corners of new Hauksson's interpolated grid
+! double precision, parameter :: UTM_X_ORIG_HAUKSSON = 122035.012d0
+! double precision, parameter :: UTM_X_END_HAUKSSON = 766968.628d0
+! double precision, parameter :: UTM_Y_ORIG_HAUKSSON = 3547232.986d0
+! double precision, parameter :: UTM_Y_END_HAUKSSON = 4098868.501d0
+
+! Lin-Shearer-Hauksson-Thurber (2007)
+! number of non-constant layers
+ integer, parameter :: NLAYERS_HAUKSSON = 8
+! depth of layers
+ double precision, parameter :: Z_HAUKSSON_LAYER_1 = 0.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_2 = -3000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_3 = -6000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_4 = -10000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_5 = -15000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_6 = -17000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_7 = -22000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_8 = -31000.d0
+
+ integer, parameter :: NGRID_NEW_HAUKSSON = 241
+
+! corners of new Hauksson's interpolated grid
+ double precision, parameter :: UTM_X_ORIG_HAUKSSON = 88021.568d0
+ double precision, parameter :: UTM_X_END_HAUKSSON = 861517.886d0
+ double precision, parameter :: UTM_Y_ORIG_HAUKSSON = 3404059.875d0
+ double precision, parameter :: UTM_Y_END_HAUKSSON = 4180234.582d0
+
+ double precision, parameter :: SPACING_UTM_X_HAUKSSON = (UTM_X_END_HAUKSSON - UTM_X_ORIG_HAUKSSON) / (NGRID_NEW_HAUKSSON-1.d0)
+ double precision, parameter :: SPACING_UTM_Y_HAUKSSON = (UTM_Y_END_HAUKSSON - UTM_Y_ORIG_HAUKSSON) / (NGRID_NEW_HAUKSSON-1.d0)
+
+! layers in the So-Cal regional model
+! DEPTH_MOHO_SOCAL = -35 km was based on Dreger and Helmberger (1990)
+! and is (July 2007) the preferred Moho depth for Dreger.
+! The depth of 32 km is used in the standard processing (Wald et al., 1995)
+! of SoCal events and is the value in the original Kanamori-Hadley (1975) model.
+ double precision, parameter :: DEPTH_5p5km_SOCAL = -5500.d0
+ double precision, parameter :: DEPTH_16km_SOCAL = -16000.d0
+ double precision, parameter :: DEPTH_MOHO_SOCAL = -32000.d0
+! double precision, parameter :: DEPTH_5p5km_SOCAL = -7700.d0
+! double precision, parameter :: DEPTH_16km_SOCAL = -18200.d0
+! double precision, parameter :: DEPTH_MOHO_SOCAL = -34200.d0
+
+! layer for Lacq (France) gas field
+ double precision, parameter :: DEPTH_INTERFACE_LACQ = -7000.d0
+
+! reference surface of the model before adding topography
+ double precision, parameter :: Z_SURFACE = 0.d0
+
+! number of points in each AVS or OpenDX quadrangular cell for movies
+ integer, parameter :: NGNOD2D_AVS_DX = 4
+
+! magic ratio for heuristic rule
+! this gives 120 degree angles in doubling
+! standard value 0.5 gives 135-135-90, which is not optimal
+ double precision, parameter :: MAGIC_RATIO = 0.6056d0
+
+! type of elements for heuristic rule
+ integer, parameter :: ITYPE_UNUSUAL_1 = 1
+ integer, parameter :: ITYPE_UNUSUAL_1p = 2
+ integer, parameter :: ITYPE_UNUSUAL_4 = 3
+ integer, parameter :: ITYPE_UNUSUAL_4p = 4
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/precision.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/precision.h (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/precision.h 2012-03-22 14:24:30 UTC (rev 19846)
@@ -0,0 +1,28 @@
+!=====================================================================
+!
+! S p e c f e m 3 D B a s i n V e r s i o n 1 . 4
+! --------------------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology July 2005
+!
+! A signed non-commercial agreement is required to use this program.
+! Please check http://www.gps.caltech.edu/research/jtromp for details.
+! Free for non-commercial academic research ONLY.
+! This program is distributed WITHOUT ANY WARRANTY whatsoever.
+! Do not redistribute this program without written permission.
+!
+!=====================================================================
+
+! precision.h. Generated from precision.h.in by configure.
+
+!
+! solver in single or double precision depending on the machine
+!
+! set to MPI_REAL to run in single precision
+! set to MPI_DOUBLE_PRECISION to run in double precision
+!
+! ALSO CHANGE FILE constants.h ACCORDINGLY
+!
+ integer, parameter :: CUSTOM_MPI_TYPE = MPI_REAL
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/smooth_sem_fun.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/smooth_sem_fun.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/smooth_sem_fun.f90 2012-03-22 14:24:30 UTC (rev 19846)
@@ -0,0 +1,317 @@
+program smooth_specfem_function
+
+! this is the embarassingly-parallel program that smooth any specfem function (primarily
+! the kernels) that has the dimension of (NGLLX,NGLLY,NGLLZ,NSPEC_MAX), notice that it
+! uses the constants.h and precision.h files from the original SEM package, and the
+! values_from_mesher.h file from the output of the mesher (or create_header_file),
+! therefore, you need to compile it for your specific case
+
+ implicit none
+ include 'mpif.h'
+ include 'constants.h'
+ include 'precision.h'
+ include 'values_from_mesher.h'
+
+! ======================================================
+
+ integer, parameter :: NSPEC_MAX=NSPEC_AB
+ integer, parameter :: NGLOB_MAX=NGLOB_AB
+
+! only include the neighboring 3 x 3 slices
+ integer, parameter :: NSLICES = 3
+ integer ,parameter :: NSLICES2 = NSLICES * NSLICES
+
+ character(len=150) :: s_nproc_xi, s_nproc_eta, s_nchunks, s_element_size, s_sigma_h, s_sigma_v
+ character(len=150) :: kernel_file_name, scratch_topo_dir, scratch_file_dir
+ integer :: sizeprocs, ier,myrank, nproc_xi, nproc_eta, nchunks, ichunk, ixi, ieta, iglob
+ integer :: islice(NSLICES2), islice0(NSLICES2), ns
+
+ real(kind=CUSTOM_REAL) :: sigma_h, sigma_h2, sigma_h3, sigma_v, sigma_v2, sigma_v3
+
+ real(kind=CUSTOM_REAL) :: x0, y0, z0, norm, norm_h, norm_v, element_size, max_old, max_new
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: factor, exp_val
+
+ character(len=150) :: ks_file, reg_name
+ character(len=150), dimension(NSLICES2) :: x_file, y_file, z_file, j_file, k_file, i_file
+
+ logical :: global_code
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_MAX) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_MAX) :: kernel, kernel_smooth
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_MAX) :: tk, bk, jacobian, xl, yl, zl, xx, yy, zz
+ real(kind=CUSTOM_REAL), dimension(NGLOB_MAX) :: x, y, z
+ real(kind=CUSTOM_REAL), dimension(NSPEC_MAX) :: cx0, cy0, cz0, cx, cy, cz
+
+ integer :: i, j, k, ispec, ii, ispec2, nspec(NSLICES2), nglob(NSLICES2)
+
+ ! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: xigll, wxgll
+ double precision, dimension(NGLLY) :: yigll, wygll
+ double precision, dimension(NGLLZ) :: zigll, wzgll
+
+ ! array with all the weights in the cube
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+
+ ! ============ program starts here =====================
+
+ ! initialize the MPI communicator and start the NPROCTOT MPI processes
+ call MPI_INIT(ier)
+ call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
+ call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
+
+ ! input arguments (nchunks=0 means basin code!)
+ call getarg(1,s_nproc_xi)
+ call getarg(2,s_nproc_eta)
+ call getarg(3,s_nchunks)
+ call getarg(4,s_element_size)
+ call getarg(5,s_sigma_h)
+ call getarg(6,s_sigma_v)
+ call getarg(7,kernel_file_name)
+ call getarg(8,scratch_file_dir)
+ call getarg(9,scratch_topo_dir)
+
+ if (trim(s_nproc_xi) == '' .or. trim(s_nproc_eta) == '' .or. trim(s_nchunks) == '' &
+ .or. trim(s_element_size) == '' .or. trim(kernel_file_name) == '' &
+ .or. trim(scratch_file_dir) == '' .or. trim(scratch_topo_dir) == '') then
+ call exit_MPI(myrank,'Usage: smooth_sem_fun nproc_xi nproc_eta nchunks element_size_on_surface(km) sigma_h(km) sigma_v(km) kernel_file_name scratch_file_dir scratch_topo_dir')
+ endif
+
+ ! read in parameter information
+ read(s_nproc_xi,*) nproc_xi
+ read(s_nproc_eta,*) nproc_eta
+ read(s_nchunks,*) nchunks
+ read(s_element_size,*) element_size
+ read(s_sigma_h,*) sigma_h
+ read(s_sigma_v,*) sigma_v
+
+ if (nchunks == 0) then
+ global_code = .false.
+ reg_name='_'
+ nchunks = 1
+ else
+ global_code = .true.
+ reg_name='_reg1_'
+ endif
+ if (sizeprocs /= nproc_xi*nproc_eta*nchunks) call exit_mpi(myrank,'Error total number of slices')
+
+ element_size = element_size * 1000 ! e.g. 9 km on the surface, 36 km at CMB
+
+ sigma_h = sigma_h * 1000.0 ! m
+ sigma_v = sigma_v * 1000.0 ! m
+
+ sigma_h2 = sigma_h ** 2
+ sigma_v2 = sigma_v ** 2
+
+ ! search radius
+ sigma_h3 = 3.0 * sigma_h + element_size
+ sigma_v3 = 3.0 * sigma_v + element_size
+
+ ! theoretic normal value
+ norm_h = 2.0*PI*sigma_h**2
+ norm_v = sqrt(2.0*PI) * sigma_v
+ norm = norm_h * norm_v
+ !norm = (sqrt(2.0*PI) * sigma) ** 3
+
+ ! GLL points
+ call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ wgll_cube(i,j,k) = wxgll(i)*wygll(j)*wzgll(k)
+ enddo
+ enddo
+ enddo
+
+! ---- figure out the neighboring 8 or 7 slices: (ichunk,ixi,ieta) index start at 0------
+ ichunk = myrank / (nproc_xi * nproc_eta)
+ ieta = (myrank - ichunk * nproc_xi * nproc_eta) / nproc_xi
+ ixi = myrank - ichunk * nproc_xi * nproc_eta - ieta * nproc_xi
+
+ ! get the neighboring slices:
+ call get_all_eight_slices(ichunk,ixi,ieta,&
+ islice0(1),islice0(2),islice0(3),islice0(4),islice0(5),islice0(6),islice0(7),islice0(8),&
+ nproc_xi,nproc_eta)
+
+ ! remove the repeated slices (only 8 for corner slices in global case)
+ islice(1) = myrank; j = 1
+ do i = 1, 8
+ if (.not. any(islice(1:i) == islice0(i)) .and. islice0(i) < sizeprocs) then
+ j = j + 1
+ islice(j) = islice0(i)
+ endif
+ enddo
+ ns = j
+
+ ! read in the topology files of the current and neighboring slices
+ do i = 1, ns
+ write(x_file(i),'(a,i6.6,a)') trim(scratch_topo_dir)//'/proc',islice(i),trim(reg_name)//'x.bin'
+ write(y_file(i),'(a,i6.6,a)') trim(scratch_topo_dir)//'/proc',islice(i),trim(reg_name)//'y.bin'
+ write(z_file(i),'(a,i6.6,a)') trim(scratch_topo_dir)//'/proc',islice(i),trim(reg_name)//'z.bin'
+ write(j_file(i),'(a,i6.6,a)') trim(scratch_topo_dir)//'/proc',islice(i),trim(reg_name)//'jacobian.bin'
+ write(i_file(i),'(a,i6.6,a)') trim(scratch_topo_dir)//'/proc',islice(i),trim(reg_name)//'ibool.bin'
+ write(k_file(i),'(a,i6.6,a)') trim(scratch_file_dir)//'/proc',islice(i),trim(reg_name)//trim(kernel_file_name)//'.bin'
+
+ nspec(i) = NSPEC_AB
+ nglob(i) = NGLOB_AB
+ enddo
+
+ ! read in myrank slice
+ open(11,file=x_file(1),status='old',form='unformatted')
+ read(11) x(1:nglob(1))
+ close(11)
+ open(11,file=y_file(1),status='old',form='unformatted')
+ read(11) y(1:nglob(1))
+ close(11)
+ open(11,file=z_file(1),status='old',form='unformatted')
+ read(11) z(1:nglob(1))
+ close(11)
+ open(11,file=i_file(1),status='old',form='unformatted')
+ read(11) ibool(:,:,:,1:nspec(1))
+ close(11)
+
+ ! get the location of the center of the elements
+ do ispec = 1, nspec(1)
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+ xl(i,j,k,ispec) = x(iglob)
+ yl(i,j,k,ispec) = y(iglob)
+ zl(i,j,k,ispec) = z(iglob)
+ enddo
+ enddo
+ enddo
+ cx0(ispec) = (xl(1,1,1,ispec) + xl(NGLLX,NGLLY,NGLLZ,ispec))/2
+ cy0(ispec) = (yl(1,1,1,ispec) + yl(NGLLX,NGLLY,NGLLZ,ispec))/2
+ cz0(ispec) = (zl(1,1,1,ispec) + zl(NGLLX,NGLLY,NGLLZ,ispec))/2
+ enddo
+
+ if (myrank == 0) write(*,*) 'start looping over elements and points for smoothing ...'
+
+ write(ks_file,'(a,i6.6,a)') trim(scratch_file_dir)//'/proc',myrank,trim(reg_name)//trim(kernel_file_name)//'_smooth.bin'
+
+ tk = 0.0; bk = 0.0; kernel_smooth=0.0
+
+ ! loop over all the slices
+ do ii = 1, ns
+
+ ! read in the topology, kernel files, calculate center of elements
+ open(11,file=x_file(ii),status='old',form='unformatted')
+ read(11) x(1:nglob(ii))
+ close(11)
+ open(11,file=y_file(ii),status='old',form='unformatted')
+ read(11) y(1:nglob(ii))
+ close(11)
+ open(11,file=z_file(ii),status='old',form='unformatted')
+ read(11) z(1:nglob(ii))
+ close(11)
+ open(11,file=i_file(ii),status='old',form='unformatted')
+ read(11) ibool(:,:,:,1:nspec(ii))
+ close(11)
+ open(11,file=j_file(ii),status='old',form='unformatted')
+ read(11) jacobian(:,:,:,1:nspec(ii))
+ close(11)
+ open(11,file=k_file(ii),status='old',form='unformatted')
+ read(11) kernel(:,:,:,1:nspec(ii))
+ close(11)
+
+! get the global maximum value of the original kernel file
+ if (ii == 1) then
+ call mpi_reduce(maxval(abs(kernel(:,:,:,1:nspec(ii)))), max_old, 1, &
+ CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+ endif
+
+ do ispec2 = 1, nspec(ii)
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec2)
+ xx(i,j,k,ispec2) = x(iglob)
+ yy(i,j,k,ispec2) = y(iglob)
+ zz(i,j,k,ispec2) = z(iglob)
+ enddo
+ enddo
+ enddo
+ cx(ispec2) = (xx(1,1,1,ispec2) + xx(NGLLX,NGLLZ,NGLLY,ispec2))/2
+ cy(ispec2) = (yy(1,1,1,ispec2) + yy(NGLLX,NGLLZ,NGLLY,ispec2))/2
+ cz(ispec2) = (zz(1,1,1,ispec2) + zz(NGLLX,NGLLZ,NGLLY,ispec2))/2
+ enddo
+
+ if (myrank == 0) write(*,*) 'slice number = ', ii
+
+ ! loop over elements to be smoothed in the current slice
+ do ispec = 1, nspec(1)
+ if (mod(ispec,100) == 0 .and. myrank == 0) write(*,*) 'ispec=', ispec
+
+ ! --- only double loop over the elements in the search radius ---
+ do ispec2 = 1, nspec(ii)
+
+ !if ( sqrt( (cx(ispec2)-cx0(ispec)) **2 + (cy(ispec2)-cy0(ispec)) ** 2 + (cz(ispec2)-cz0(ispec)) ** 2) > sigma3) cycle
+ if ( sqrt( (cx(ispec2)-cx0(ispec)) **2 + (cy(ispec2)-cy0(ispec)) ** 2 ) > sigma_h3 .or. sqrt( (cz(ispec2)-cz0(ispec)) ** 2) > sigma_v3 ) cycle
+
+ factor(:,:,:) = jacobian(:,:,:,ispec2) * wgll_cube(:,:,:) ! integration factors
+
+ ! loop over GLL points of the elements in current slice (ispec)
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ x0 = xl(i,j,k,ispec); y0 = yl(i,j,k,ispec); z0 = zl(i,j,k,ispec) ! current point (i,j,k,ispec)
+
+ !exp_val(:,:,:) = exp( -((xx(:,:,:,ispec2)-x0)**2+(yy(:,:,:,ispec2)-y0)**2 &
+ ! +(zz(:,:,:,ispec2)-z0)**2 )/(2*sigma2) )*factor(:,:,:)
+
+ exp_val(:,:,:) = exp( -(xx(:,:,:,ispec2)-x0)**2/(2.0*sigma_h2) &
+ -(yy(:,:,:,ispec2)-y0)**2/(2.0*sigma_h2) &
+ -(zz(:,:,:,ispec2)-z0)**2/(2.0*sigma_v2) ) * factor(:,:,:)
+
+ tk(i,j,k,ispec) = tk(i,j,k,ispec) + sum(exp_val(:,:,:) * kernel(:,:,:,ispec2))
+ bk(i,j,k,ispec) = bk(i,j,k,ispec) + sum(exp_val(:,:,:))
+
+ enddo
+ enddo
+ enddo ! (i,j,k)
+ enddo ! (ispec2)
+ enddo ! (ispec)
+ enddo ! islice
+
+ if (myrank == 0) write(*,*) 'Done with integration ...'
+
+ ! compute the smoothed kernel values
+ do ispec = 1, nspec(1)
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ if (bk(i,j,k,ispec) < 0.01 * norm) then ! check the normalization criterion
+ print *, 'Problem here --- ', myrank, ispec, i, j, k, bk(i,j,k,ispec), 0.01 * norm
+ call exit_mpi(myrank, 'Error computing Gaussian function on the grid')
+ endif
+
+ kernel_smooth(i,j,k,ispec) = tk(i,j,k,ispec)/bk(i,j,k,ispec)
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ open(11,file=trim(ks_file),status='unknown',form='unformatted')
+ ! Note: output the following instead of kernel_smooth(:,:,:,1:nspec(1)) to create files of the same sizes
+ write(11) kernel_smooth(:,:,:,:)
+ close(11)
+
+ ! the maximum value for the smoothed kernel
+ call mpi_reduce(maxval(abs(kernel_smooth(:,:,:,1:nspec(1)))), max_new, 1, &
+ CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+
+ if (myrank == 0) then
+ print *, 'Maximum data value before smoothing = ', max_old
+ print *, 'Maximum data value after smoothing = ', max_new
+ endif
+
+! stop all the MPI processes, and exit
+ call MPI_FINALIZE(ier)
+
+end program smooth_specfem_function
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/values_from_mesher.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/values_from_mesher.h (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/smooth/src/values_from_mesher.h 2012-03-22 14:24:30 UTC (rev 19846)
@@ -0,0 +1,54 @@
+
+ !
+ ! this is the parameter file for static compilation of the solver
+ !
+ ! mesh statistics:
+ ! ---------------
+ !
+ ! number of processors = 168
+ !
+ ! number of ES nodes = 21.00000
+ ! percentage of total 640 ES nodes = 3.281250 %
+ ! total memory available on these ES nodes (Gb) = 336.0000
+ !
+ ! max points per processor = max vector length = 163941
+ ! min vector length = 25
+ ! min critical vector length = 75
+ !
+ ! on ES and SX-5, make sure "loopcnt=" parameter
+ ! in Makefile is greater than 163941
+ !
+ ! total elements per AB slice = 2412
+ ! total points per AB slice = 163941
+ !
+ ! total for full mesh:
+ ! -------------------
+ !
+ ! exact total number of spectral elements in entire mesh =
+ ! 405216
+ ! approximate total number of points in entire mesh =
+ ! 27542088.0000000
+ ! approximate total number of degrees of freedom in entire mesh =
+ ! 82626264.0000000
+ !
+ ! resolution of the mesh at the surface:
+ ! -------------------------------------
+ !
+ ! spectral elements along X = 336
+ ! spectral elements along Y = 288
+ ! GLL points along X = 1345
+ ! GLL points along Y = 1153
+ ! average distance between points along X in m = 475.4175
+ ! average distance between points along Y in m = 436.8476
+ !
+
+ integer, parameter :: NSPEC_AB = 2412
+ integer, parameter :: NGLOB_AB = 163941
+
+ !
+ ! number of time steps = 27300
+ !
+ integer, parameter :: NSPEC_ATTENUATION = 1
+ logical, parameter :: ATTENUATION_VAL = .false.
+ logical, parameter :: ANISOTROPY_VAL = .false.
+
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/Makefile
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/Makefile 2012-03-22 04:17:38 UTC (rev 19845)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/cluster/sum_kernel/src/Makefile 2012-03-22 14:24:30 UTC (rev 19846)
@@ -3,25 +3,15 @@
OBJS = sum_kernels.o exit_mpi.o
-OBJS_G = sum_kernels_globe.o exit_mpi.o
+all: sum_kernels
-OBJS_G_PRE = sum_preconditioned_kernels_globe.o exit_mpi.o
-
-all: sum_kernels sum_kernels_globe sum_preconditioned_kernels_globe
-
sum_kernels: $(OBJS)
$(FC) $(FFLAGS) -o $@ $(OBJS)
-sum_kernels_globe: $(OBJS_G)
- $(FC) $(FFLAGS) -o $@ $(OBJS_G)
-
-sum_preconditioned_kernels_globe: $(OBJS_G_PRE)
- $(FC) $(FFLAGS) -o $@ $(OBJS_G_PRE)
-
.SUFFIXES: .o .f90
.f90.o:
$(FC) $(FFLAGS) -c $<
clean:
- rm -rf *.o *~ sum_kernels_globe sum_preconditioned_kernels_globe sum_kernels
+ rm -rf *.o *~ sum_kernels
More information about the CIG-COMMITS
mailing list