[cig-commits] r20154 - in seismo/3D/SPECFEM3D/trunk: . src/generate_databases src/shared src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Wed May 16 12:08:48 PDT 2012
Author: danielpeter
Date: 2012-05-16 12:08:47 -0700 (Wed, 16 May 2012)
New Revision: 20154
Added:
seismo/3D/SPECFEM3D/trunk/src/shared/sum_kernels.f90
Modified:
seismo/3D/SPECFEM3D/trunk/Makefile.in
seismo/3D/SPECFEM3D/trunk/src/generate_databases/model_external_values.f90
seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in
Log:
adds sum_kernels.f90 for kernel summations; compile with: make xsum_kernels
Modified: seismo/3D/SPECFEM3D/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/Makefile.in 2012-05-16 18:17:20 UTC (rev 20153)
+++ seismo/3D/SPECFEM3D/trunk/Makefile.in 2012-05-16 19:08:47 UTC (rev 20154)
@@ -68,10 +68,6 @@
@COND_PYRE_FALSE@ specfem3D \
@COND_PYRE_FALSE@ $(EMPTY_MACRO)
-# @COND_PYRE_FALSE@ combine_vol_data \
-# @COND_PYRE_FALSE@ combine_surf_data \
-# @COND_PYRE_FALSE@ convolve_source_timefunction \
-# @COND_PYRE_FALSE@ smooth_vol_data \
default: $(DEFAULT)
@@ -98,6 +94,7 @@
combine_vol_data: xcombine_vol_data
combine_surf_data: xcombine_surf_data
smooth_vol_data: xsmooth_vol_data
+sum_kernels: xsum_kernels
model_update: xmodel_update
check_mesh_quality_CUBIT_Abaqus: xcheck_mesh_quality_CUBIT_Abaqus
@@ -152,6 +149,9 @@
xsmooth_vol_data: required reqspec
(cd src/specfem3D ; make xsmooth_vol_data)
+xsum_kernels: required reqspec
+ (cd src/specfem3D ; make xsum_kernels)
+
xmodel_update: required reqspec xspecfem3D
(cd src/specfem3D ; make xmodel_update)
@@ -177,6 +177,7 @@
@echo " xcombine_vol_data"
@echo " xcombine_surf_data"
@echo " xsmooth_vol_data"
+ @echo " xsum_kernels"
@echo " xmodel_update"
@echo " xcheck_mesh_quality_CUBIT_Abaqus"
@echo ""
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/model_external_values.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/model_external_values.f90 2012-05-16 18:17:20 UTC (rev 20153)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/model_external_values.f90 2012-05-16 19:08:47 UTC (rev 20154)
@@ -156,6 +156,9 @@
y = ymesh
z = zmesh
+ ! note: z coordinate will be negative below surface
+ ! convention is z-axis points up
+
! model dimensions
xmin = 0. ! minval(xstore_dummy)
xmax = 134000. ! maxval(xstore_dummy)
@@ -168,21 +171,22 @@
call get_topo_elevation_free_closest(x,y,elevation,distmin, &
nspec,nglob_dummy,ibool,xstore_dummy,ystore_dummy,zstore_dummy, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
+
! depth in Z-direction
if( distmin < HUGEVAL ) then
depth = elevation - z
else
- depth = zmax - z
+ depth = zmin - z
endif
! normalizes depth between 0 and 1
if( abs( zmax - zmin ) > TINYVAL ) depth = depth / (zmax - zmin)
- ! super-imposes values
- !rho = 2.6910d0+0.6924d0*depth
- !vp = 4.1875d0+3.9382d0*depth
- !vs = 2.1519d0+2.3481d0*depth
+ ! initial values (in m/s and kg/m^3)
+ rho = 2691.0d0
+ vp = 4187.5d0
+ vs = 2151.9d0
! adds a velocity depth gradient
! (e.g. from PREM mantle gradients:
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90 2012-05-16 18:17:20 UTC (rev 20153)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90 2012-05-16 19:08:47 UTC (rev 20154)
@@ -24,6 +24,27 @@
!
!=====================================================================
+ module vtk
+ !-------------------------------------------------------------
+ ! USER PARAMETER
+
+ ! outputs as VTK ASCII file
+ logical,parameter :: USE_VTK_OUTPUT = .true.
+
+ !-------------------------------------------------------------
+
+ ! global point data
+ real,dimension(:),allocatable :: total_dat
+
+ ! maximum number of slices
+ integer,parameter :: MAX_NUM_NODES = 600
+
+ end module vtk
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
program combine_paraview_data_ext_mesh
! puts the output of SPECFEM3D into '***.mesh' format,
@@ -36,6 +57,7 @@
!
! works for external, unregular meshes
+ use vtk
implicit none
include 'constants.h'
@@ -51,8 +73,12 @@
integer :: NSPEC_AB, NGLOB_AB
integer :: numpoin
+
integer :: i, ios, it, ier
- integer :: iproc, proc1, proc2, num_node, node_list(2000)
+ integer :: iproc, proc1, proc2, num_node
+
+ integer,dimension(MAX_NUM_NODES) :: node_list
+
integer :: np, ne, npp, nee, nelement, njunk
character(len=256) :: sline, arg(6), filename, indir, outdir
@@ -76,7 +102,7 @@
logical :: ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
character(len=256) LOCAL_PATH
integer :: IMODEL
-
+
! checks given arguments
print *
print *,'Recombining ParaView data for slices'
@@ -117,6 +143,7 @@
read(sline,*,iostat=ios) njunk
if (ios /= 0) exit
num_node = num_node + 1
+ if( num_node > MAX_NUM_NODES ) stop 'error number of slices exceeds MAX_NUM_NODES...'
node_list(num_node) = njunk
enddo
close(20)
@@ -128,7 +155,9 @@
read(arg(1),*) proc1
read(arg(2),*) proc2
do iproc = proc1, proc2
- node_list(iproc - proc1 + 1) = iproc
+ it = iproc - proc1 + 1
+ if( it > MAX_NUM_NODES ) stop 'error number of slices exceeds MAX_NUM_NODES...'
+ node_list(it) = iproc
enddo
num_node = proc2 - proc1 + 1
filename = arg(3)
@@ -157,10 +186,21 @@
print *, 'Slice list: '
print *, node_list(1:num_node)
- ! open paraview output mesh file
- mesh_file = trim(outdir) // '/' // trim(filename)//'.mesh'
- call open_file(trim(mesh_file)//char(0))
+ if( USE_VTK_OUTPUT ) then
+ mesh_file = trim(outdir) // '/' // trim(filename)//'.vtk'
+ open(IOVTK,file=mesh_file(1:len_trim(mesh_file)),status='unknown',iostat=ios)
+ if( ios /= 0 ) stop 'error opening vtk output file'
+ write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOVTK,'(a)') 'material model VTK file'
+ write(IOVTK,'(a)') 'ASCII'
+ write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ else
+ ! open paraview output mesh file
+ mesh_file = trim(outdir) // '/' // trim(filename)//'.mesh'
+ call open_file(trim(mesh_file)//char(0))
+ endif
+
! counts total number of points (all slices)
npp = 0
nee = 0
@@ -229,11 +269,11 @@
if (.not. HIGH_RESOLUTION_MESH) then
! writes out element corners only
call cvd_write_corners(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore,dat, &
- it,npp,numpoin)
+ it,npp,numpoin,np)
else
! high resolution, all GLL points
call cvd_write_GLL_points(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore,dat,&
- it,npp,numpoin)
+ it,npp,numpoin,np)
endif
print*,' points:',np,numpoin
@@ -246,6 +286,8 @@
enddo ! all slices for points
+ if( USE_VTK_OUTPUT) write(IOVTK,*) ""
+
if (np /= npp) stop 'Error: Number of total points are not consistent'
print *, 'Total number of points: ', np
print *, ' '
@@ -293,6 +335,8 @@
enddo ! num_node
+ if( USE_VTK_OUTPUT) write(IOVTK,*) ""
+
! checks with total number of elements
if (ne /= nee) then
print*,'error: number of elements counted:',ne,'total:',nee
@@ -300,9 +344,25 @@
endif
print *, 'Total number of elements: ', ne
- ! close mesh file
- call close_file()
+ if( USE_VTK_OUTPUT) then
+ ! type: hexahedrons
+ write(IOVTK,'(a,i12)') "CELL_TYPES ",nee
+ write(IOVTK,*) (12,it=1,nee)
+ write(IOVTK,*) ""
+ write(IOVTK,'(a,i12)') "POINT_DATA ",npp
+ write(IOVTK,'(a)') "SCALARS "//trim(filename)//" float"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do it = 1,npp
+ write(IOVTK,*) total_dat(it)
+ enddo
+ write(IOVTK,*) ""
+ close(IOVTK)
+ else
+ ! close mesh file
+ call close_file()
+ endif
+
print *, 'Done writing '//trim(mesh_file)
end program combine_paraview_data_ext_mesh
@@ -317,11 +377,14 @@
! counts total number of points and elements for external meshes in given slice list
! returns: total number of elements (nee) and number of points (npp)
+ use vtk
implicit none
include 'constants.h'
- integer,intent(in) :: num_node,node_list(300)
+ integer,intent(in) :: num_node
+ integer,dimension(MAX_NUM_NODES),intent(in) :: node_list
character(len=256),intent(in) :: LOCAL_PATH
+
integer,intent(out) :: npp,nee
logical,intent(in) :: HIGH_RESOLUTION_MESH
@@ -416,10 +479,10 @@
subroutine cvd_write_corners(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore,dat,&
- it,npp,numpoin)
+ it,npp,numpoin,np)
! writes out locations of spectral element corners only
-
+ use vtk
implicit none
include 'constants.h'
@@ -428,7 +491,7 @@
real(kind=CUSTOM_REAL),dimension(NGLOB_AB) :: xstore, ystore, zstore
real,dimension(NGLLY,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: dat
integer:: it
- integer :: npp,numpoin
+ integer :: npp,numpoin,np
! local parameters
logical,dimension(:),allocatable :: mask_ibool
@@ -438,7 +501,15 @@
! writes out total number of points
if (it == 1) then
- call write_integer(npp)
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK, '(a,i12,a)') 'POINTS ', npp, ' float'
+ ! creates array to hold point data
+ allocate(total_dat(npp),stat=ier)
+ if( ier /= 0 ) stop 'error allocating total dat array'
+ total_dat(:) = 0.0
+ else
+ call write_integer(npp)
+ endif
endif
! writes our corner point locations
@@ -461,21 +532,31 @@
x = xstore(iglob1)
y = ystore(iglob1)
z = zstore(iglob1)
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(1,1,1,ispec))
- mask_ibool(iglob1) = .true.
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK,'(3e18.6)') x,y,z
+ total_dat(np+numpoin) = dat(1,1,1,ispec)
+ else
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(1,1,1,ispec))
+ endif
+ mask_ibool(iglob1) = .true.
endif
if(.not. mask_ibool(iglob2)) then
numpoin = numpoin + 1
x = xstore(iglob2)
y = ystore(iglob2)
z = zstore(iglob2)
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(NGLLX,1,1,ispec))
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK,'(3e18.6)') x,y,z
+ total_dat(np+numpoin) = dat(NGLLX,1,1,ispec)
+ else
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(NGLLX,1,1,ispec))
+ endif
mask_ibool(iglob2) = .true.
endif
if(.not. mask_ibool(iglob3)) then
@@ -483,10 +564,15 @@
x = xstore(iglob3)
y = ystore(iglob3)
z = zstore(iglob3)
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(NGLLX,NGLLY,1,ispec))
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK,'(3e18.6)') x,y,z
+ total_dat(np+numpoin) = dat(NGLLX,NGLLY,1,ispec)
+ else
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(NGLLX,NGLLY,1,ispec))
+ endif
mask_ibool(iglob3) = .true.
endif
if(.not. mask_ibool(iglob4)) then
@@ -494,10 +580,15 @@
x = xstore(iglob4)
y = ystore(iglob4)
z = zstore(iglob4)
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(1,NGLLY,1,ispec))
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK,'(3e18.6)') x,y,z
+ total_dat(np+numpoin) = dat(1,NGLLY,1,ispec)
+ else
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(1,NGLLY,1,ispec))
+ endif
mask_ibool(iglob4) = .true.
endif
if(.not. mask_ibool(iglob5)) then
@@ -505,10 +596,15 @@
x = xstore(iglob5)
y = ystore(iglob5)
z = zstore(iglob5)
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(1,1,NGLLZ,ispec))
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK,'(3e18.6)') x,y,z
+ total_dat(np+numpoin) = dat(1,1,NGLLZ,ispec)
+ else
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(1,1,NGLLZ,ispec))
+ endif
mask_ibool(iglob5) = .true.
endif
if(.not. mask_ibool(iglob6)) then
@@ -516,10 +612,15 @@
x = xstore(iglob6)
y = ystore(iglob6)
z = zstore(iglob6)
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(NGLLX,1,NGLLZ,ispec))
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK,'(3e18.6)') x,y,z
+ total_dat(np+numpoin) = dat(NGLLX,1,NGLLZ,ispec)
+ else
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(NGLLX,1,NGLLZ,ispec))
+ endif
mask_ibool(iglob6) = .true.
endif
if(.not. mask_ibool(iglob7)) then
@@ -527,10 +628,15 @@
x = xstore(iglob7)
y = ystore(iglob7)
z = zstore(iglob7)
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(NGLLX,NGLLY,NGLLZ,ispec))
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK,'(3e18.6)') x,y,z
+ total_dat(np+numpoin) = dat(NGLLX,NGLLY,NGLLZ,ispec)
+ else
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(NGLLX,NGLLY,NGLLZ,ispec))
+ endif
mask_ibool(iglob7) = .true.
endif
if(.not. mask_ibool(iglob8)) then
@@ -538,10 +644,15 @@
x = xstore(iglob8)
y = ystore(iglob8)
z = zstore(iglob8)
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(1,NGLLY,NGLLZ,ispec))
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK,'(3e18.6)') x,y,z
+ total_dat(np+numpoin) = dat(1,NGLLY,NGLLZ,ispec)
+ else
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(1,NGLLY,NGLLZ,ispec))
+ endif
mask_ibool(iglob8) = .true.
endif
enddo ! ispec
@@ -553,10 +664,10 @@
subroutine cvd_write_GLL_points(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore,dat,&
- it,npp,numpoin)
+ it,npp,numpoin,np)
! writes out locations of all GLL points of spectral elements
-
+ use vtk
implicit none
include 'constants.h'
@@ -564,7 +675,7 @@
integer,dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: ibool
real(kind=CUSTOM_REAL),dimension(NGLOB_AB) :: xstore, ystore, zstore
real,dimension(NGLLY,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: dat
- integer:: it,npp,numpoin
+ integer:: it,npp,numpoin,np
! local parameters
logical,dimension(:),allocatable :: mask_ibool
@@ -573,7 +684,15 @@
! writes out total number of points
if (it == 1) then
- call write_integer(npp)
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK, '(a,i12,a)') 'POINTS ', npp, ' float'
+ ! creates array to hold point data
+ allocate(total_dat(npp),stat=ier)
+ if( ier /= 0 ) stop 'error allocating total dat array'
+ total_dat(:) = 0.0
+ else
+ call write_integer(npp)
+ endif
endif
! writes out point locations and values
@@ -592,10 +711,15 @@
x = xstore(iglob)
y = ystore(iglob)
z = zstore(iglob)
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(i,j,k,ispec))
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK,'(3e18.6)') x,y,z
+ total_dat(np+numpoin) = dat(i,j,k,ispec)
+ else
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(i,j,k,ispec))
+ endif
mask_ibool(iglob) = .true.
endif
enddo ! i
@@ -612,6 +736,7 @@
subroutine cvd_write_corner_elements(NSPEC_AB,NGLOB_AB,ibool,&
np,nelement,it,nee,numpoin)
+ use vtk
implicit none
include 'constants.h'
@@ -628,7 +753,12 @@
! outputs total number of elements for all slices
if (it == 1) then
- call write_integer(nee)
+ if( USE_VTK_OUTPUT ) then
+ ! note: indices for vtk start at 0
+ write(IOVTK,'(a,i12,i12)') "CELLS ",nee,nee*9
+ else
+ call write_integer(nee)
+ endif
end if
! writes out element indices
@@ -702,14 +832,18 @@
n7 = num_ibool(iglob7) -1 + np
n8 = num_ibool(iglob8) -1 + np
- call write_integer(n1)
- call write_integer(n2)
- call write_integer(n3)
- call write_integer(n4)
- call write_integer(n5)
- call write_integer(n6)
- call write_integer(n7)
- call write_integer(n8)
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK,'(9i12)') 8,n1,n2,n3,n4,n5,n6,n7,n8
+ else
+ call write_integer(n1)
+ call write_integer(n2)
+ call write_integer(n3)
+ call write_integer(n4)
+ call write_integer(n5)
+ call write_integer(n6)
+ call write_integer(n7)
+ call write_integer(n8)
+ endif
enddo
@@ -729,7 +863,7 @@
np,nelement,it,nee,numpoin)
! writes out indices of elements given by GLL points
-
+ use vtk
implicit none
include 'constants.h'
@@ -746,8 +880,13 @@
! outputs total number of elements for all slices
if (it == 1) then
- !nee = nelement * num_node
- call write_integer(nee)
+ if( USE_VTK_OUTPUT ) then
+ ! note: indices for vtk start at 0
+ write(IOVTK,'(a,i12,i12)') "CELLS ",nee,nee*9
+ else
+ !nee = nelement * num_node
+ call write_integer(nee)
+ endif
endif
! sets numbering num_ibool respecting mask
@@ -794,14 +933,19 @@
n6 = num_ibool(iglob6)+np-1
n7 = num_ibool(iglob7)+np-1
n8 = num_ibool(iglob8)+np-1
- call write_integer(n1)
- call write_integer(n2)
- call write_integer(n3)
- call write_integer(n4)
- call write_integer(n5)
- call write_integer(n6)
- call write_integer(n7)
- call write_integer(n8)
+
+ if( USE_VTK_OUTPUT ) then
+ write(IOVTK,'(9i12)') 8,n1,n2,n3,n4,n5,n6,n7,n8
+ else
+ call write_integer(n1)
+ call write_integer(n2)
+ call write_integer(n3)
+ call write_integer(n4)
+ call write_integer(n5)
+ call write_integer(n6)
+ call write_integer(n7)
+ call write_integer(n8)
+ endif
enddo
enddo
enddo
Added: seismo/3D/SPECFEM3D/trunk/src/shared/sum_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/sum_kernels.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/sum_kernels.f90 2012-05-16 19:08:47 UTC (rev 20154)
@@ -0,0 +1,461 @@
+! sum_preconditioned_kernels
+!
+! this program can be used for event kernel summation,
+! where it sums up transverse isotropic kernel files:
+!
+! - proc***_reg1_bulk_c_kernel.bin
+! - proc***_reg1_bulk_betav_kernel.bin
+! - proc***_reg1_bulk_betah_kernel.bin
+! - proc***_reg1_eta_kernel.bin
+!
+! this version uses the approximate Hessian stored as
+! - proc***_reg1_hess_kernel.bin
+! to precondition the transverse isotropic kernel files
+!
+! input file: kernels_list.txt
+! lists all event kernel directories which should be summed together
+!
+! input directory: INPUT_KERNELS/
+! contains links to all event kernel directories (listed in "kernels_list.txt")
+!
+! output directory: OUTPUT_SUM/
+! the resulting kernel files will be stored in this directory
+
+module sum_par
+
+ include 'constants.h'
+
+ ! USER PARAMETERS
+
+ ! by default, this algorithm uses transverse isotropic (bulk,bulk_betav,bulk_betah,eta) kernels to sum up
+ ! if you prefer using isotropic kernels, set flags below accordingly
+
+ ! if you prefer using isotropic kernels (bulk,bulk_beta,rho) kernels, set this flag to true
+ logical, parameter :: USE_ISO_KERNELS = .false.
+
+ ! if you prefer isotropic (alpha,beta,rho) kernels, set this flag to true
+ logical, parameter :: USE_ALPHA_BETA_RHO = .false.
+
+ ! 1 permille of maximum for inverting hessian
+ real(kind=CUSTOM_REAL),parameter :: THRESHOLD_HESS = 1.e-3
+
+ ! sums all hessians before inverting and preconditioning
+ logical, parameter :: USE_HESS_SUM = .true.
+
+ ! uses source mask to blend out source elements
+ logical, parameter :: USE_SOURCE_MASK = .false.
+
+ ! maximum number of kernels listed
+ integer, parameter :: MAX_NUM_NODES = 1000
+
+ ! default list name
+ character(len=150),parameter :: kernel_file_list = 'kernels_list.txt'
+
+ ! mesh size
+ integer :: NSPEC_AB, NGLOB_AB
+
+end module sum_par
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+program sum_kernels
+
+ use sum_par
+ implicit none
+
+ include 'mpif.h'
+ include 'precision.h'
+
+
+ character(len=150) :: kernel_list(MAX_NUM_NODES)
+ character(len=150) :: sline, kernel_name,prname_lp
+ integer :: nker, myrank, sizeprocs, ier
+ integer :: ios
+
+ ! for read_parameter_files
+ double precision :: DT
+ double precision :: HDUR_MOVIE
+ integer :: NPROC,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP, &
+ UTM_PROJECTION_ZONE,SIMULATION_TYPE
+ integer :: NSOURCES,NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY
+ integer :: NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
+ integer :: IMODEL
+ logical :: MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
+ USE_HIGHRES_FOR_MOVIES,SUPPRESS_UTM_PROJECTION
+ logical :: ATTENUATION,USE_OLSEN_ATTENUATION, &
+ OCEANS,TOPOGRAPHY
+ logical :: ABSORBING_CONDITIONS,SAVE_FORWARD
+ logical :: ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
+ character(len=256) LOCAL_PATH
+
+ ! ============ 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)
+
+ if(myrank==0) then
+ write(*,*) 'sum_preconditioned_kernels:'
+ write(*,*)
+ write(*,*) 'reading kernel list: '
+ endif
+ call mpi_barrier(MPI_COMM_WORLD,ier)
+
+ ! reads in event list
+ nker=0
+ open(unit = 20, file = trim(kernel_file_list), status = 'old',iostat = ios)
+ if (ios /= 0) then
+ print *,'Error opening ',trim(kernel_file_list),myrank
+ stop 1
+ endif
+ do while (1 == 1)
+ read(20,'(a)',iostat=ios) sline
+ if (ios /= 0) exit
+ nker = nker+1
+ kernel_list(nker) = sline
+ enddo
+ close(20)
+ if( myrank == 0 ) then
+ write(*,*) ' ',nker,' events'
+ write(*,*)
+ endif
+
+ ! needs local_path for mesh files
+ call read_parameter_file( NPROC,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,DT, &
+ UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION, &
+ ATTENUATION,USE_OLSEN_ATTENUATION,LOCAL_PATH,NSOURCES, &
+ OCEANS,TOPOGRAPHY,ANISOTROPY,ABSORBING_CONDITIONS, &
+ MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
+ NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
+ SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION, &
+ NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
+ NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY,IMODEL)
+
+ ! checks if number of MPI process as specified
+ if (sizeprocs /= NPROC) then
+ if( myrank == 0 ) then
+ print*,''
+ print*,'error: run xsum_kernels with the same number of MPI processes '
+ print*,' as specified in Par_file by NPROC when slices were created'
+ print*,''
+ print*,'for example: mpirun -np ',NPROC,' ./xsum_kernels ...'
+ print*,''
+ endif
+ call exit_mpi(myrank,'Error total number of slices')
+ endif
+ call mpi_barrier(MPI_COMM_WORLD,ier)
+
+ ! reads mesh file
+ !
+ ! needs to get array dimensions
+
+ ! opens external mesh file
+ write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'/proc',myrank,'_'//'external_mesh.bin'
+ open(unit=27,file=trim(prname_lp),&
+ status='old',action='read',form='unformatted',iostat=ios)
+ if( ier /= 0 ) then
+ print*,'error: could not open database '
+ print*,'path: ',trim(prname_lp)
+ call exit_mpi(myrank, 'error reading external mesh file')
+ endif
+
+ ! gets number of elements and global points for this partition
+ read(27) NSPEC_AB
+ read(27) NGLOB_AB
+
+ close(27)
+
+ ! user output
+ if(myrank == 0) then
+ print*
+ print*,' rank:',myrank,' summing kernels'
+ print*,kernel_list(1:nker)
+ print*
+ endif
+
+ ! synchronizes
+ call mpi_barrier(MPI_COMM_WORLD,ier)
+
+ ! sums up kernels
+ if( USE_ISO_KERNELS ) then
+
+ if( myrank == 0 ) write(*,*) 'isotropic kernels: bulk_c, bulk_beta, rho'
+
+ ! isotropic kernels
+ kernel_name = 'reg1_bulk_c_kernel'
+ call sum_kernel_pre(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_bulk_beta_kernel'
+ call sum_kernel_pre(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_rho_kernel'
+ call sum_kernel_pre(kernel_name,kernel_list,nker,myrank)
+
+ else if( USE_ALPHA_BETA_RHO ) then
+
+ ! isotropic kernels
+ if( myrank == 0 ) write(*,*) 'isotropic kernels: alpha, beta, rho'
+
+ kernel_name = 'reg1_alpha_kernel'
+ call sum_kernel_pre(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_beta_kernel'
+ call sum_kernel_pre(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_rho_kernel'
+ call sum_kernel_pre(kernel_name,kernel_list,nker,myrank)
+
+ else
+
+ ! transverse isotropic kernels
+ if( myrank == 0 ) write(*,*) 'transverse isotropic kernels: bulk_c, bulk_betav, bulk_betah,eta'
+
+ kernel_name = 'reg1_bulk_c_kernel'
+ call sum_kernel_pre(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_bulk_betav_kernel'
+ call sum_kernel_pre(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_bulk_betah_kernel'
+ call sum_kernel_pre(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_eta_kernel'
+ call sum_kernel_pre(kernel_name,kernel_list,nker,myrank)
+
+ endif
+
+ if(myrank==0) write(*,*) 'done writing all kernels'
+
+ ! stop all the MPI processes, and exit
+ call MPI_FINALIZE(ier)
+
+end program sum_kernels
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine sum_kernel_pre(kernel_name,kernel_list,nker,myrank)
+
+ use sum_par
+ implicit none
+
+ include 'mpif.h'
+ include 'precision.h'
+
+ real(kind=CUSTOM_REAL) :: norm,norm_sum
+ character(len=150) :: kernel_name,kernel_list(MAX_NUM_NODES)
+ integer :: nker,myrank,ier
+
+ ! local parameters
+ character(len=150) :: k_file
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: &
+ kernel,hess,total_kernel
+ integer :: iker,ios
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: total_hess,mask_source
+
+ ! initializes arrays
+ allocate(kernel(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ hess(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ total_kernel(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+
+
+ if( USE_HESS_SUM ) then
+ allocate( total_hess(NGLLX,NGLLY,NGLLZ,NSPEC_AB) )
+ total_hess(:,:,:,:) = 0.0_CUSTOM_REAL
+ endif
+
+ if( USE_SOURCE_MASK ) then
+ allocate( mask_source(NGLLX,NGLLY,NGLLZ,NSPEC_AB) )
+ mask_source(:,:,:,:) = 1.0_CUSTOM_REAL
+ endif
+
+ ! loops over all event kernels
+ total_kernel = 0._CUSTOM_REAL
+ do iker = 1, nker
+ ! user output
+ if(myrank==0) then
+ write(*,*) 'reading in event kernel for: ',trim(kernel_name)
+ write(*,*) ' and preconditioner: ','reg1_hess_kernel'
+ write(*,*) ' ',iker, ' out of ', nker
+ endif
+
+ ! sensitivity kernel / frechet derivative
+ kernel = 0._CUSTOM_REAL
+ write(k_file,'(a,i6.6,a)') 'INPUT_KERNELS/'//trim(kernel_list(iker)) &
+ //'/proc',myrank,'_'//trim(kernel_name)//'.bin'
+
+ open(12,file=trim(k_file),status='old',form='unformatted',action='read',iostat=ios)
+ if( ios /= 0 ) then
+ write(*,*) ' kernel not found:',trim(k_file)
+ cycle
+ endif
+ read(12) kernel
+ close(12)
+
+ ! outputs norm of kernel
+ norm = sum( kernel * kernel )
+ call mpi_reduce(norm,norm_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ print*,' norm kernel: ',sqrt(norm_sum)
+ endif
+
+
+ ! approximate Hessian
+ hess = 0._CUSTOM_REAL
+ write(k_file,'(a,i6.6,a)') 'INPUT_KERNELS/'//trim(kernel_list(iker)) &
+ //'/proc',myrank,'_reg1_hess_kernel.bin'
+
+ open(12,file=trim(k_file),status='old',form='unformatted',action='read',iostat=ios)
+ if( ios /= 0 ) then
+ write(*,*) ' hess not found:',trim(k_file)
+ else
+ read(12) hess
+ close(12)
+ endif
+
+ ! outputs norm of preconditioner
+ norm = sum( hess * hess )
+ call mpi_reduce(norm,norm_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ print*,' norm preconditioner: ',sqrt(norm_sum)
+ print*
+ endif
+
+ ! note: we take absolute values for hessian (as proposed by Yang)
+ hess = abs(hess)
+
+ ! source mask
+ if( USE_SOURCE_MASK ) then
+ ! reads in mask
+ write(k_file,'(a,i6.6,a)') 'INPUT_KERNELS/'//trim(kernel_list(iker)) &
+ //'/proc',myrank,'_reg1_mask_source.bin'
+ open(12,file=trim(k_file),status='old',form='unformatted',action='read',iostat=ios)
+ if( ios /= 0 ) then
+ write(*,*) ' file not found:',trim(k_file)
+ cycle
+ endif
+ read(12) mask_source
+ close(12)
+
+ ! masks source elements
+ kernel = kernel * mask_source
+
+ endif
+
+ ! precondition
+ if( USE_HESS_SUM ) then
+
+ ! sums up hessians first
+ total_hess = total_hess + hess
+
+ else
+
+ ! inverts hessian
+ call invert_hess( myrank,hess )
+
+ ! preconditions each event kernel with its hessian
+ kernel = kernel * hess
+
+ endif
+
+ ! sums all kernels from each event
+ total_kernel = total_kernel + kernel
+
+ enddo
+
+ ! preconditions summed kernels with summed hessians
+ if( USE_HESS_SUM ) then
+
+ ! inverts hessian matrix
+ call invert_hess( myrank,total_hess )
+
+ ! preconditions kernel
+ total_kernel = total_kernel * total_hess
+
+ endif
+
+ ! stores summed kernels
+ if(myrank==0) write(*,*) 'writing out summed kernel for: ',trim(kernel_name)
+
+ write(k_file,'(a,i6.6,a)') 'OUTPUT_SUM/proc',myrank,'_'//trim(kernel_name)//'.bin'
+
+ open(12,file=trim(k_file),form='unformatted',status='unknown',action='write',iostat=ios)
+ if( ios /= 0 ) then
+ write(*,*) 'ERROR kernel not written:',trim(k_file)
+ stop 'error kernel write'
+ endif
+ write(12) total_kernel
+ close(12)
+
+ if(myrank==0) write(*,*)
+
+ ! frees memory
+ deallocate(kernel,hess,total_kernel)
+ if( USE_HESS_SUM ) deallocate(total_hess)
+ if( USE_SOURCE_MASK ) deallocate(mask_source)
+
+end subroutine sum_kernel_pre
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine invert_hess( myrank,hess_matrix )
+
+! inverts the hessian matrix
+! the approximate hessian is only defined for diagonal elements: like
+! H_nn = \frac{ \partial^2 \chi }{ \partial \rho_n \partial \rho_n }
+! on all GLL points, which are indexed (i,j,k,ispec)
+
+ use sum_par
+ implicit none
+
+ include 'mpif.h'
+ include 'precision.h'
+
+ integer :: myrank
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ hess_matrix
+
+ ! local parameters
+ real(kind=CUSTOM_REAL) :: maxh,maxh_all
+ integer :: ier
+
+ ! maximum value of hessian
+ maxh = maxval( abs(hess_matrix) )
+
+ ! determines maximum from all slices on master
+ call mpi_allreduce(maxh,maxh_all,1,CUSTOM_MPI_TYPE,MPI_MAX,MPI_COMM_WORLD,ier)
+
+ ! user output
+ if( myrank == 0 ) then
+ print*
+ print*,'hessian maximum: ',maxh_all
+ print*
+ endif
+
+ ! normalizes hessian
+ if( maxh_all < 1.e-18 ) then
+ ! hessian is zero, re-initializes
+ hess_matrix = 1.0_CUSTOM_REAL
+ !call exit_mpi(myrank,'error hessian too small')
+ else
+ ! since hessian has absolute values, this scales between [0,1]
+ hess_matrix = hess_matrix / maxh_all
+ endif
+
+
+ ! inverts hessian values
+ where( abs(hess_matrix(:,:,:,:)) > THRESHOLD_HESS )
+ hess_matrix = 1.0_CUSTOM_REAL / hess_matrix
+ elsewhere
+ hess_matrix = 1.0_CUSTOM_REAL / THRESHOLD_HESS
+ endwhere
+
+ ! rescales hessian
+ !hess_matrix = hess_matrix * maxh_all
+
+end subroutine invert_hess
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in 2012-05-16 18:17:20 UTC (rev 20153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in 2012-05-16 19:08:47 UTC (rev 20154)
@@ -162,7 +162,6 @@
$O/get_value_parameters.o \
$O/param_reader.o \
$O/exit_mpi.o \
- $O/parallel.o \
$O/read_mesh_databases.o \
$O/create_name_database.o \
$O/check_mesh_resolution.o \
@@ -172,6 +171,17 @@
$O/write_VTK_data.o \
$(EMPTY_MACRO)
+
+SUM_KERNELS_OBJECTS = \
+ $O/sum_kernels.o \
+ $O/read_parameter_file.o \
+ $O/read_value_parameters.o \
+ $O/get_value_parameters.o \
+ $O/param_reader.o \
+ $O/exit_mpi.o \
+ $(EMPTY_MACRO)
+
+
# objects toggled between the parallel and serial version
@COND_MPI_TRUE at COND_MPI_OBJECTS = $O/parallel.o
@COND_MPI_FALSE at COND_MPI_OBJECTS = $O/serial.o
@@ -194,6 +204,7 @@
@COND_PYRE_FALSE@ combine_surf_data \
@COND_PYRE_FALSE@ convolve_source_timefunction \
@COND_PYRE_FALSE@ smooth_vol_data \
+ at COND_PYRE_FALSE@ sum_kernels \
@COND_PYRE_FALSE@ model_update \
@COND_PYRE_FALSE@ $(EMPTY_MACRO)
@@ -222,6 +233,7 @@
combine_vol_data: xcombine_vol_data
combine_surf_data: xcombine_surf_data
smooth_vol_data: xsmooth_vol_data
+sum_kernels: xsum_kernels
model_update: xmodel_update
xconvolve_source_timefunction: $O/convolve_source_timefunction.o
@@ -239,14 +251,18 @@
xsmooth_vol_data: $O/smooth_vol_data.o $O/read_parameter_file.o $O/read_value_parameters.o $O/get_value_parameters.o $O/param_reader.o $O/gll_library.o $O/exit_mpi.o $O/parallel.o
${FCLINK} -o ${E}/xsmooth_vol_data $O/smooth_vol_data.o $O/read_parameter_file.o $O/read_value_parameters.o $O/get_value_parameters.o $O/param_reader.o $O/gll_library.o $O/exit_mpi.o $O/parallel.o $(MPILIBS)
-xmodel_update: $(MODEL_UPD_OBJECTS)
- ${FCLINK} -o ${E}/xmodel_update $(MODEL_UPD_OBJECTS) $(MPILIBS)
+xsum_kernels: $(SUM_KERNELS_OBJECTS) $(COND_MPI_OBJECTS)
+ ${FCLINK} -o ${E}/xsum_kernels $(SUM_KERNELS_OBJECTS) $(COND_MPI_OBJECTS) $(MPILIBS)
+xmodel_update: $(MODEL_UPD_OBJECTS) $(COND_MPI_OBJECTS)
+ ${FCLINK} -o ${E}/xmodel_update $(MODEL_UPD_OBJECTS) $(COND_MPI_OBJECTS) $(MPILIBS)
+
clean:
rm -f $O/* *.o *.gnu *.mod $(OUTPUT)/timestamp* $(OUTPUT)/starttime*txt work.pc* \
xgenerate_databases xspecfem3D \
xconvolve_source_timefunction \
- xcreate_movie_shakemap_AVS_DX_GMT xcombine_vol_data xcombine_surf_data xsmooth_vol_data xmodel_update
+ xcreate_movie_shakemap_AVS_DX_GMT xcombine_vol_data xcombine_surf_data \
+ xsmooth_vol_data xmodel_update xsum_kernels
###
### rule for the archive library
@@ -547,6 +563,13 @@
${MPIFCCOMPILE_NO_CHECK} -c -o $O/smooth_vol_data.o ${SHARED}/smooth_vol_data.f90
##
+## kernel summation
+##
+
+$O/sum_kernels.o: $(SHARED)constants.h ${SHARED}/sum_kernels.f90
+ ${MPIFCCOMPILE_NO_CHECK} -c -o $O/sum_kernels.o ${SHARED}/sum_kernels.f90
+
+##
## model update
##
More information about the CIG-COMMITS
mailing list