adds a tomography tool xsum_kernels; updates routine exit_MPI() to be…

… similar to 3D version

EXAMPLES/Tape2007_kernel/adj_seismogram_Tape2007.f90

@@ -64,6 +64,7 @@ program adj_seismogram
   double precision :: seism_veloc(NSTEP),seism_accel(NSTEP),ft_bar(NSTEP)
   character(len=3) :: compr(2),comp(3)
   character(len=150) :: filename
+  integer :: ier
 
   NDIM=3
   comp = (/"BXX","BXY","BXZ"/)
@@ -87,8 +88,9 @@ program adj_seismogram
 
      do icomp = 1, NDIMr
 
-        filename = 'OUTPUT_FILES/'//trim(station_name(irec))//'.AA.'// compr(icomp) // '.semd'
-        open(unit = 10, file = trim(filename))
+        filename = 'OUTPUT_FILES/'//'AA.'//trim(station_name(irec))//'.'// compr(icomp) // '.semd'
+        open(unit = 10, file = trim(filename),status='old',iostat=ier)
+        if (ier /= 0) stop 'Error opening trace file'
 
         do itime = 1,NSTEP
            read(10,*) time , seism(itime,icomp)
@@ -119,8 +121,9 @@ program adj_seismogram
 
         print *,comp(icomp)
 
-        filename = 'SEM/'//trim(station_name(irec))//'.AA.'// comp(icomp) // '.adj'
-        open(unit = 11, file = trim(filename))
+        filename = 'SEM/'//'AA.'//trim(station_name(irec))//'.'// comp(icomp) // '.adj'
+        open(unit = 11, file = trim(filename),status='unknown',iostat=ier)
+        if (ier /= 0) stop 'Error opening SEM adjoint source file'
 
         time_window(:) = 0.d0
         seism_win(:) = seism(:,icomp)
@@ -170,7 +173,7 @@ program adj_seismogram
 
   enddo
   print *,'*************************'
-  print *,'The input files (S****.AA.BXX/BXY/BXZ.adj) needed to run the adjoint simulation are in SEM'
+  print *,'The input files (AA.S****.BXX/BXY/BXZ.adj) needed to run the adjoint simulation are in SEM'
   print *,'*************************'
 
 end program adj_seismogram

Makefile.in

@@ -211,6 +211,7 @@ SUBDIRS = \
 	shared \
 	specfem2D \
 	tomography/postprocess_sensitivity_kernels \
+	tomography \
 	$(EMPTY_MACRO)
 
 # default targets for the pure Fortran version
@@ -225,7 +226,7 @@ DEFAULT = \
 
 default: $(DEFAULT)
 
-all: default auxiliaries postprocess
+all: default auxiliaries postprocess tomography
 
 
 ifdef CLEAN
@@ -267,6 +268,9 @@ help:
 	@echo "    xcombine_sem"
 	@echo "    xsmooth_sem"
 	@echo ""
+	@echo "- tomography tools: [make tomography]"
+	@echo "    xsum_kernels"
+	@echo ""
 
 
 .PHONY: all default backup clean realclean help

configure

@@ -7435,7 +7435,7 @@ $as_echo "## ----------------------------------- ##
 ## ----------------------------------- ##"
 
 # Output results.
-ac_config_files="$ac_config_files Makefile setup/constants.h setup/precision.h"
+ac_config_files="$ac_config_files Makefile setup/constants.h setup/constants_tomography.h setup/precision.h"
 
 
 if test "$USE_BUNDLED_SCOTCH" = 1; then
@@ -8176,6 +8176,7 @@ do
     "setup/config.h") CONFIG_HEADERS="$CONFIG_HEADERS setup/config.h" ;;
     "Makefile") CONFIG_FILES="$CONFIG_FILES Makefile" ;;
     "setup/constants.h") CONFIG_FILES="$CONFIG_FILES setup/constants.h" ;;
+    "setup/constants_tomography.h") CONFIG_FILES="$CONFIG_FILES setup/constants_tomography.h" ;;
     "setup/precision.h") CONFIG_FILES="$CONFIG_FILES setup/precision.h" ;;
     "${SCOTCH_DIR}/src/Makefile.inc") CONFIG_FILES="$CONFIG_FILES ${SCOTCH_DIR}/src/Makefile.inc" ;;
     "bin") CONFIG_COMMANDS="$CONFIG_COMMANDS bin" ;;

configure.ac

@@ -320,6 +320,7 @@ AS_BOX([setting up default simulation setup])
 AC_CONFIG_FILES([
 	Makefile
 	setup/constants.h
+	setup/constants_tomography.h
 	setup/precision.h
 ])
 

setup/constants.h.in

@@ -235,6 +235,9 @@
 ! do not use tags for MPI messages, use dummy tag instead
   integer, parameter :: itag = 0
 
+! for the Gauss-Lobatto-Legendre points and weights
+  double precision, parameter :: GAUSSALPHA = 0.d0,GAUSSBETA = 0.d0
+
 ! number of lines per source in SOURCE file
   integer, parameter :: NLINES_PER_SOURCE = 14
 
@@ -246,6 +249,9 @@
 ! ignore variable name field (junk) at the beginning of each input line
   logical, parameter :: IGNORE_JUNK = .true., DONT_IGNORE_JUNK = .false.
 
+! maximum length of strings used for paths, reading from files, etc.
+  integer, parameter :: MAX_STRING_LEN = 512
+
 !!-----------------------------------------------------------
 !!
 !! DOMAINS

setup/constants_tomography.h.in

@@ -0,0 +1,131 @@
+!=====================================================================
+!
+!                   S P E C F E M 2 D  Version 7 . 0
+!
+!=====================================================================
+
+!
+!--- user can modify parameters below
+!
+
+!!-----------------------------------------------------------
+!!
+!! Model parameterization
+!!
+!!-----------------------------------------------------------
+! by default, this algorithm uses isotropic (alpha,beta) 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 = .true.
+
+  ! region code
+  ! (local version "_", globe version "_reg1_" for crust/mantle
+  character(len=8),parameter :: REG = "_"
+
+!!-----------------------------------------------------------
+!!
+!! Directory structure
+!!
+!!-----------------------------------------------------------
+
+  ! directory containing kernels from current iteration
+  character(len=MAX_STRING_LEN) :: INPUT_KERNELS_DIR = 'INPUT_GRADIENT/'
+
+  ! directory containing kernels from previous iteration
+  character(len=MAX_STRING_LEN) :: KERNEL_OLD_DIR = './KERNELS/OUTPUT_SUM.old'
+
+  ! directory containing model database files
+  character(len=MAX_STRING_LEN) :: INPUT_MODEL_DIR = 'INPUT_MODEL/'
+
+  ! directory containing external_mesh.bin database files
+  character(len=MAX_STRING_LEN) :: INPUT_DATABASES_DIR = 'topo/'
+
+  ! destination for newly updated model database files
+  character(len=MAX_STRING_LEN) :: OUTPUT_MODEL_DIR = 'OUTPUT_MODEL/'
+
+  ! destination for log files containing minimum, maximum vales
+  character(len=MAX_STRING_LEN) :: OUTPUT_STATISTICS_DIR = 'OUTPUT_MODEL/'
+
+  ! file containing absolute or relative paths to kernel directories
+  character(len=MAX_STRING_LEN) :: KERNEL_FILE_LIST = './kernels_list.txt'
+
+  ! produce logs containing minimum, maximum vales and other statistics
+  logical :: PRINT_STATISTICS_FILES = .false.
+
+  ! maximum number of kernel paths
+  integer, parameter :: MAX_KERNEL_PATHS = 10000
+
+!!-----------------------------------------------------------
+!!
+!! Scaling laws
+!!
+!!-----------------------------------------------------------
+  ! ignore rho kernel, but use density perturbations as a scaling of Vs perturbations
+  logical, parameter :: USE_RHO_SCALING = .true.
+
+  ! in case of rho scaling, specifies density scaling factor with shear perturbations
+  ! see e.g. Montagner & Anderson (1989), Panning & Romanowicz (2006)
+  real(kind=CUSTOM_REAL),parameter :: RHO_SCALING = 0.33_CUSTOM_REAL
+
+!!-----------------------------------------------------------
+!!
+!! Transversely isotropic (TI) model constraints
+!!
+!!-----------------------------------------------------------
+  ! constraint on eta model
+  real(kind=CUSTOM_REAL),parameter :: LIMIT_ETA_MIN = 0.5_CUSTOM_REAL
+  real(kind=CUSTOM_REAL),parameter :: LIMIT_ETA_MAX = 1.5_CUSTOM_REAL
+
+!!-----------------------------------------------------------
+!!
+!! Approximate hessian
+!!
+!!-----------------------------------------------------------
+  ! 1 permille of maximum for inverting hessian
+  real(kind=CUSTOM_REAL),parameter :: THRESHOLD_HESS = 1.e-3
+
+  ! sums all hessians before inverting and preconditioning
+  ! by default should be set to .true.
+  logical, parameter :: USE_HESS_SUM = .true.
+
+!!-----------------------------------------------------------
+!!
+!! Maximum kernel scaling
+!!
+!!-----------------------------------------------------------
+! kernel values are maximum at very shallow depth (due to receivers at surface) which leads to strong
+! model updates closest to the surface. scaling the kernel values, such that the maximum is taken slightly below
+! the surface (between 1km - 5km) leads to a "more balanced" gradient, i.e., a better model update in deeper parts
+
+  ! by default, sets maximum update in this depth range
+  logical,parameter :: USE_DEPTH_RANGE_MAXIMUM = .false.
+
+  ! depths
+  ! top at 1km depth
+  real(kind=CUSTOM_REAL),parameter :: R_TOP = -1000.0 ! shallow depth
+  ! bottom at 5km depth
+  real(kind=CUSTOM_REAL),parameter :: R_BOTTOM = -5000.0 ! deep depth
+
+!!-----------------------------------------------------------
+!!
+!! Source mask
+!!
+!!-----------------------------------------------------------
+  ! uses source mask to blend out source elements
+  logical, parameter :: USE_SOURCE_MASK = .false.
+
+
+!!-----------------------------------------------------------
+!!
+!! Deprecated
+!!
+!!-----------------------------------------------------------
+
+  ! determines separate step length for each material parameter in nonlinear conjugate gradient model update
+  logical,parameter :: USE_SEPARATE_CG_STEPLENGTHS = .false.
+
+

src/auxiliaries/rules.mk

@@ -87,6 +87,7 @@ check_quality_external_mesh_auxiliaries_OBJECTS = \
 check_quality_external_mesh_auxiliaries_SHARED_OBJECTS = \
 	$O/define_shape_functions.shared.o \
 	$O/exit_mpi.shared.o \
+	$O/parallel.shared.o \
 	$O/param_reader.cc.o \
 	$O/read_value_parameters.shared.o \
 	$(EMPTY_MACRO)