Input/Output files
TomoATT requires 3 input files :
- input parameter file (setup file for simulation parameters)
- source receiver file (source and receiver definitions and observation arrival times)
- initial model file (3d initial model)
It is recommended to use TomoATT's python API for creating the template project directory and input files. This library is available from here (opens in a new tab).
If PyPI is not available on your environment, you may git clone the repository and install it manually by,
git clone --branch devel https://github.com/MIGG-NTU/PyTomoATT/tree/devel
cd PyTomoATT
pip install -e .After instaling PyTomoATT, you can create a template project directory with the following command
pta init_pjt path/to/project_dirBelow is the explanation of each input file.
1. input parameter file
All the necessary parameters for setuping a calculation are described in input parameter file in yaml format (opens in a new tab).
Below is an example of input parameter file for making a forward simulation.
version: 3
#################################################
# computational domian #
#################################################
domain:
min_max_dep: [-10, 10] # depth in km
min_max_lat: [37.7, 42.3] # latitude in degree
min_max_lon: [22.7, 27.3] # longitude in degree
n_rtp: [10, 50, 50] # number of nodes in depth,latitude,longitude direction
#################################################
# traveltime data file path #
#################################################
source:
src_rec_file: OUTPUT_FILES/src_rec_file_forward.dat # source receiver file path
swap_src_rec: true # swap source and receiver
#################################################
# initial model file path #
#################################################
model:
init_model_path: ./test_model_init.h5 # path to initial model file
model_1d_name: dummy_model_1d_name # 1D model name used in teleseismic 2D solver (iasp91, ak135, user_defined is available), defined in include/1d_model.h
#################################################
# parallel computation settings #
#################################################
parallel: # parameters for parallel computation
n_sims: 1 # number of simultanoues runs (parallel the sources)
ndiv_rtp: [1, 1, 1] # number of subdivision on each direction (parallel the computional domain)
nproc_sub: 2 # number of processors for sweep parallelization (parallel the fast sweep method)
use_gpu: false # true if use gpu (EXPERIMENTAL)
############################################
# output file setting #
############################################
output_setting:
output_dir: ./OUTPUT_FILES/ # path to output director (default is ./OUTPUT_FILES/)
output_source_field: true # output the calculated field of all sources
output_model_dat: false # output model_parameters_inv_0000.dat (data in text format) or not.
output_final_model: true # output merged final model (final_model.h5) or not.
output_in_process: true # output model at each inv iteration or not.
output_in_process_data: true # output src_rec_file at each inv iteration or not.
single_precision_output: false # output results in single precision or not.
verbose_output_level: 0 # output internal parameters, if no, only model parameters are out.
output_file_format: 0
#################################################
# inversion or forward modeling #
#################################################
# run mode
# 0 for forward simulation only,
# 1 for inversion
# 2 for earthquake relocation
# 3 for inversion + earthquake relocation
run_mode: 1
###################################################
# model update parameters setting #
###################################################
model_update:
max_iterations: 3 # maximum number of inversion iterations
optim_method: 0 # optimization method. 0 : grad_descent, 1 : halve-stepping, 2 : lbfgs (EXPERIMENTAL)
#common parameters for all optim methods
step_length: 0.01 # step length of model perturbation at each iteration. 0.01 means maximum 1% perturbation for each iteration.
# parameters for optim_method 0 (gradient_descent)
optim_method_0:
step_length_decay: 0.9 # if objective function increase, step size -> step length * step_length_decay. default: 0.9
# parameters for optim_method 1 (halve-stepping) or 2 (lbfgs)
optim_method_1_2:
max_sub_iterations: 20 # maximum number of each sub-iteration
regularization_weight: 0.5 # weight value for regularization (lbfgs mode only)
coefs_regulalization_rtp: [1, 1, 1] # regularization coefficients for rtp (lbfgs mode only)
# smoothing
smoothing:
smooth_method: 0 # 0: multiparametrization, 1: laplacian smoothing (EXPERIMENTAL)
l_smooth_rtp: [1, 0.0174533, 0.0174533] # smoothing coefficients for laplacian smoothing
# parameters for smooth method 0 (multigrid model parametrization)
# inversion grid can be viewed in OUTPUT_FILES/inversion_grid.txt
n_inversion_grid: 5 # number of inversion grid sets
# inversion grid type
type_invgrid_dep: 0 # 0: uniform inversion grid, 1: flexible grid
type_invgrid_lat: 0 # 0: uniform inversion grid, 1: flexible grid
type_invgrid_lon: 0 # 0: uniform inversion grid, 1: flexible grid
# settings for uniform inversion grid (if type_*_inv : 0)
n_inv_dep_lat_lon: [5, 10, 10] # number of the base inversion grid points
min_max_dep_inv: [-10, 10] # depth in km (Radius of the earth is defined in config.h/R_earth)
min_max_lat_inv: [37.7, 42.3] # latitude in degree
min_max_lon_inv: [22.7, 27.3] # longitude in degree
# settings for flexible inversion grid (if type_*_inv : 1)
dep_inv: [1, 1, 1]
lat_inv: [1, 1, 1]
lon_inv: [1, 1, 1]
# if we want to use another inversion grid for inverting anisotropy, set invgrid_ani: true (default: false)
invgrid_ani: false
# settings for flexible inversion grid for anisotropy (only flexible grid input is provided)
# anisotropic inversion grid type (used only if invgrid_ani : true)
type_invgrid_dep_ani: 0 # 0: uniform inversion grid, 1: flexible grid
type_invgrid_lat_ani: 0 # 0: uniform inversion grid, 1: flexible grid
type_invgrid_lon_ani: 0 # 0: uniform inversion grid, 1: flexible grid
# settings for uniform anisotropic inversion grid (used only if invgrid_ani : true)
n_inv_dep_lat_lon_ani: [1, 1, 1] # number of the base inversion grid points
min_max_dep_inv_ani: [-99999, -99999] # depth in km (Radius of the earth is defined in config.h/R_earth)
min_max_lat_inv_ani: [-99999, -99999] # latitude in degree
min_max_lon_inv_ani: [-99999, -99999] # longitude in degree
# settings for flexible inversion grid for anisotropy (used only if invgrid_ani : true)
dep_inv_ani: [1, 1, 1]
lat_inv_ani: [1, 1, 1]
lon_inv_ani: [1, 1, 1]
# inversion grid volume rescale (kernel -> kernel / volume of inversion grid mesh),
# this precondition may be carefully applied if the sizes of inversion grids are unbalanced
invgrid_volume_rescale: true
# path to station correction file (under development)
use_sta_correction: false
# sta_correction_file: dummy_sta_correction_file # station correction file path
step_length_sc: 0.001 # step length relate to the update of station correction terms
# In the following data subsection, XXX_weight means a weight is assigned to the data, influencing the objective function and gradient
# XXX_weight : [d1,d2,w1,w2] means:
# if XXX < d1, weight = w1
# if d1 <= XXX < d2, weight = w1 + (XXX-d1)/(d2-d1)*(w2-w1), (linear interpolation)
# if d2 <= XXX , weight = w2
# You can easily set w1 = w2 = 1.0 to normalize the weight related to XXX.
# -------------- using absolute traveltime data --------------
abs_time:
use_abs_time: true # 'true' for using absolute traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
residual_weight: [1, 3, 1, 1] # XXX is the absolute traveltime residual (second) = abs(t^{obs}_{n,i} - t^{syn}_{n,j})
distance_weight: [50, 150, 1, 1] # XXX is epicenter distance (km) between the source and receiver related to the data
# -------------- using common source differential traveltime data --------------
cs_dif_time:
use_cs_time: false # 'true' for using common source differential traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
residual_weight: [1, 3, 1, 0.1] # XXX is the common source differential traveltime residual (second) = abs(t^{obs}_{n,i} - t^{obs}_{n,j} - t^{syn}_{n,i} + t^{syn}_{n,j}).
azimuthal_weight: [15, 30, 1, 0.1] # XXX is the azimuth difference between two separate stations related to the common source.
# -------------- using common receiver differential traveltime data --------------
cr_dif_time:
use_cr_time: false # 'true' for using common receiver differential traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
residual_weight: [1, 3, 1, 0.1] # XXX is the common receiver differential traveltime residual (second) = abs(t^{obs}_{n,i} - t^{obs}_{m,i} - t^{syn}_{n,i} + t^{syn}_{m,i})
azimuthal_weight: [15, 30, 1, 0.1] # XXX is the azimuth difference between two separate sources related to the common receiver.
# -------------- global weight of different types of data (to balance the weight of different data) --------------
global_weight:
balance_data_weight: true # true: over the total weight of the each type of the data. false: use original weight (below weight for each type of data needs to be set)
abs_time_weight: 1 # weight of absolute traveltime data after balance, default: 1.0
cs_dif_time_local_weight: 1 # weight of common source differential traveltime data after balance, default: 1.0
cr_dif_time_local_weight: 1 # weight of common receiver differential traveltime data after balance, default: 1.0
teleseismic_weight: 1 # weight of teleseismic data after balance, default: 1.0 (exclude in this version)
# -------------- inversion parameters --------------
update_slowness : true # update slowness (velocity) or not. default: true
update_azi_ani : true # update azimuthal anisotropy (xi, eta) or not. default: false
# -------------- for teleseismic inversion (under development) --------------
# depth_taper : [d1,d2] means:
# if XXX < d1, kernel <- kernel * 0.0
# if d1 <= XXX < d2, kernel <- kernel * (XXX-d1)/(d2-d1), (linear interpolation)
# if d2 <= XXX , kernel <- kernel * 1.0
# You can easily set d1 = -200, d1 = -100 to remove this taper.
depth_taper : [-200, -100]
#################################################
# relocation parameters setting #
#################################################
relocation: # update earthquake hypocenter and origin time (when run_mode : 2 and 3)
min_Ndata: 4 # if the number of data of the earthquake is less than <min_Ndata>, the earthquake will not be relocated. defaut value: 4
# relocation_strategy
step_length : 0.01 # step length of relocation perturbation at each iteration. 0.01 means maximum 1% perturbation for each iteration.
step_length_decay : 0.9 # if objective function increase, step size -> step length * step_length_decay. default: 0.9
rescaling_dep_lat_lon_ortime : [10, 10, 10, 1] # The perturbation is related to <rescaling_dep_lat_lon_ortime>. Unit: km,km,km,second
max_change_dep_lat_lon_ortime : [5, 5, 5, 0.5] # the change of dep,lat,lon,ortime do not exceed max_change. Unit: km,km,km,second
max_iterations : 100 # maximum number of iterations for relocation
tol_gradient : 0.0001 # if the norm of gradient is smaller than the tolerance, the iteration of relocation terminates
# -------------- using absolute traveltime data --------------
abs_time:
use_abs_time : true # 'true' for using absolute traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
residual_weight : [1, 3, 1, 0.1] # XXX is the absolute traveltime residual (second) = abs(t^{obs}_{n,i} - t^{syn}_{n,j})
distance_weight : [50, 150, 1, 0.1] # XXX is epicenter distance (km) between the source and receiver related to the data
# -------------- using common receiver differential traveltime data --------------
cr_dif_time:
use_cr_time : true # 'true' for using common receiver differential traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
residual_weight : [1, 3, 1, 0.1] # XXX is the common receiver differential traveltime residual (second) = abs(t^{obs}_{n,i} - t^{obs}_{m,i} - t^{syn}_{n,i} + t^{syn}_{m,i})
azimuthal_weight : [10, 30, 1, 0.1] # XXX is the azimuth difference between two separate sources related to the common receiver.
# -------------- global weight of different types of data (to balance the weight of different data) --------------
global_weight:
balance_data_weight: false # true: over the total weight of the each type of the data. no: use original weight (below weight for each type of data needs to be set)
abs_time_local_weight: 1 # weight of absolute traveltime data for relocation after balance, default: 1.0
cr_dif_time_local_weight: 1 # weight of common receiver differential traveltime data for relocation after balance, default: 1.0
####################################################################
# inversion strategy for tomography and relocation #
####################################################################
inversion_strategy: # update model parameters and earthquake hypocenter iteratively (when run_mode : 3)
inv_mode : 0 # 0 for update model parameters and relocation iteratively. (other options for future work)
# for inv_mode : 0, parameters below are required
inv_mode_0: # update model for <model_update_N_iter> steps, then update location for <relocation_N_iter> steps, and repeat the process for <max_loop> loops.
model_update_N_iter : 1
relocation_N_iter : 1
max_loop : 10
####################################################################
# Eikonal solver setups #
####################################################################
calculation:
convergence_tolerance: 0.0001 # threshold value for checking the convergence for each forward/adjoint run
max_iterations: 500 # number of maximum iteration for each forward/adjoint run
stencil_order: 3 # order of stencil, 1 or 3
stencil_type: 0 # 0: , 1: first-order upwind scheme (only sweep_type 0 is supported)
sweep_type: 1 # 0: legacy, 1: cuthill-mckee with shm parallelization
There are categories and sub-categories with setup parameters. Below is the explanation of each category and sub-category. If a category tag or sub-category tag is missing from the input parameter file, the default value will be used.
Computational domain
Settings for a global domain size and number of grid points of the simulation.
domain:
min_max_dep: [-10, 10] # depth in km
min_max_lat: [37.7, 42.3] # latitude in degree
min_max_lon: [22.7, 27.3] # longitude in degree
n_rtp: [10, 50, 50] # number of nodes in depth,latitude,longitude directionTraveltime data file path
Specification of input source receiver file and swap source and receiver.
A file path to source receiver file is set with src_rec_file tag.
If swap_src_rec is set to true, source and receiver will be swapped, which is useful when the number of receiver is smaller than the number of sources.
source:
src_rec_file: OUTPUT_FILES/src_rec_file_forward.dat # source receiver file path
swap_src_rec: true # swap source and receiveswap_src_rec parameter is ignored for teleseimic sources.
Inital model settings
Input model file path is used for the calculation.
model_1d_name is used for 2D teleseismic solver. iasp91 and ak135 are predefined 1D models.
user_defined is available for user defined 1D model, which may be defined in include/1d_models.h/model_1d_prem,
which has 2 column format, depth (km) and velocity (km/s).
model:
init_model_path: ./test_model_init.h5 # path to initial model file
model_1d_name: dummy_model_1d_name # 1D model name used in teleseismic 2D solver (iasp91, ak135, user_defined is available), defined in include/1d_model.hParallel computation settings
The number of simultaneous runs are set with n_sims tag. Total number of sources are divided into n_sims groups, and each group is calculated simultaneously.
ndiv_rtp is the number of subdivision on each direction.
nproc_sub is the number of processors for sweep parallelization.
parallel: # parameters for parallel computation
n_sims: 1 # number of simultanoues runs (parallel the sources)
ndiv_rtp: [1, 1, 1] # number of subdivision on each direction (parallel the computional domain)
nproc_sub: 2 # number of processors for sweep parallelization (parallel the fast sweep method)
use_gpu: false # true if use gpu (EXPERIMENTAL)The RAM space is allocated for each simultaneous run group, while the full source receiver data is stored only by the first processor of the first simultaneous run group. The other run groups retain only the sources which are calculated within the group.
Sweeping plane parallelization (nproc_sub >= 2) is only available for stencil_type is 0, otherwise this parameter needs to be set to 1.
Output file setting
Output directory path is set with output_dir tag, which stores all the output files.
If output_source_field is true, calculated field e.g. traveltime field, adjoint field will be output for all the sources, otherwise only the model parameters will be output.
output_model_dat == true enforce the output of model parameters in text format, which is useful for debugging.
output_final_model == true output the final model in separated file. No post-processing (treatment of the ghost layers) is required to merge the model parameters.
If output_in_process is true, the calculated model parameters at each iteration will be output.
If output_in_process_data is true, the source receiver file at each iteration will be output.
single_precision_output == true output the results in single precision, which is useful for reducing the output file size.
verbose_output_level == 0 output only the model parameters, == 1 outputs some additional parameters e.g. the kernels before/after smoothing etc.
output_file_format == 0 output the results in HDF5, == 1 output the results in ASCII format.
output_setting:
output_dir: ./OUTPUT_FILES/ # path to output director (default is ./OUTPUT_FILES/)
output_source_field: true # output the calculated field of all sources
output_model_dat: false # output model_parameters_inv_0000.dat (data in text format) or not.
output_final_model: true # output merged final model (final_model.h5) or not.
output_in_process: true # output model at each inv iteration or not.
output_in_process_data: true # output src_rec_file at each inv iteration or not.
single_precision_output: false # output results in single precision or not.
verbose_output_level: 0 # output internal parameters, if no, only model parameters are out.
output_file_format: 0 # 0: HDF5, 1: ASCIIRun mode
# run mode
# 0 for forward simulation only,
# 1 for inversion
# 2 for earthquake relocation
# 3 for inversion + earthquake relocation
run_mode: 1Model update parameters setting
This section is for selecting the optimization method, setting the maximum number of iterations, and setting the step length of model perturbation.
model_update:
max_iterations: 3 # maximum number of inversion iterations
optim_method: 0 # optimization method. 0 : grad_descent, 1 : halve-stepping, 2 : lbfgs (EXPERIMENTAL)
#common parameters for all optim methods
step_length: 0.01 # step length of model perturbation at each iteration. 0.01 means maximum 1% perturbation for each iteration.
# parameters for optim_method 0 (gradient_descent)
optim_method_0:
step_length_decay: 0.9 # if objective function increase, step size -> step length * step_length_decay. default: 0.9
# parameters for optim_method 1 (halve-stepping) or 2 (lbfgs)
optim_method_1_2:
max_sub_iterations: 20 # maximum number of each sub-iteration
regularization_weight: 0.5 # weight value for regularization (lbfgs mode only)
coefs_regulalization_rtp: [1, 1, 1] # regularization coefficients for rtp (lbfgs mode only)
# smoothing
smoothing:
smooth_method: 0 # 0: multiparametrization, 1: laplacian smoothing (EXPERIMENTAL)
l_smooth_rtp: [1, 0.0174533, 0.0174533] # smoothing coefficients for laplacian smoothingThe settings below is for configuring the inversion grid for smooth_method == 0.
Equal spacing/flexible spacing grid type is selected with type_invgrid_dep, type_invgrid_lat, type_invgrid_lon.
# parameters for smooth method 0 (multigrid model parametrization)
# inversion grid can be viewed in OUTPUT_FILES/inversion_grid.txt
n_inversion_grid: 5 # number of inversion grid sets
# inversion grid type
type_invgrid_dep: 0 # 0: uniform inversion grid, 1: flexible grid
type_invgrid_lat: 0 # 0: uniform inversion grid, 1: flexible grid
type_invgrid_lon: 0 # 0: uniform inversion grid, 1: flexible grid
# settings for uniform inversion grid (if type_*_inv : 0)
n_inv_dep_lat_lon: [5, 10, 10] # number of the base inversion grid points
min_max_dep_inv: [-10, 10] # depth in km (Radius of the earth is defined in config.h/R_earth)
min_max_lat_inv: [37.7, 42.3] # latitude in degree
min_max_lon_inv: [22.7, 27.3] # longitude in degree
# settings for flexible inversion grid (if type_*_inv : 1)
dep_inv: [1, 1, 1]
lat_inv: [1, 1, 1]
lon_inv: [1, 1, 1]By invgrid_ani == true, the different inversion grid may be defined for anisotropic parameters.
# if we want to use another inversion grid for inverting anisotropy, set invgrid_ani: true (default: false)
invgrid_ani: false
# settings for flexible inversion grid for anisotropy (only flexible grid input is provided)
# anisotropic inversion grid type (used only if invgrid_ani : true)
type_invgrid_dep_ani: 0 # 0: uniform inversion grid, 1: flexible grid
type_invgrid_lat_ani: 0 # 0: uniform inversion grid, 1: flexible grid
type_invgrid_lon_ani: 0 # 0: uniform inversion grid, 1: flexible grid
# settings for uniform anisotropic inversion grid (used only if invgrid_ani : true)
n_inv_dep_lat_lon_ani: [1, 1, 1] # number of the base inversion grid points
min_max_dep_inv_ani: [-99999, -99999] # depth in km (Radius of the earth is defined in config.h/R_earth)
min_max_lat_inv_ani: [-99999, -99999] # latitude in degree
min_max_lon_inv_ani: [-99999, -99999] # longitude in degree
# settings for flexible inversion grid for anisotropy (used only if invgrid_ani : true)
dep_inv_ani: [1, 1, 1]
lat_inv_ani: [1, 1, 1]
lon_inv_ani: [1, 1, 1]
# inversion grid volume rescale (kernel -> kernel / volume of inversion grid mesh),
# this precondition may be carefully applied if the sizes of inversion grids are unbalanced
invgrid_volume_rescale: true
# path to station correction file (under development)
use_sta_correction: false
# sta_correction_file: dummy_sta_correction_file # station correction file path
step_length_sc: 0.001 # step length relate to the update of station correction termsThis section is for specifying how the traveltime data in the source receiver file is used during the inversion.
The linear weight on residual (difference between synthetic and observed traveltime) is set with residual_weight tag.
The linear weight on distance (between source and receiver) may also be set distance_weight.
# In the following data subsection, XXX_weight means a weight is assigned to the data, influencing the objective function and gradient
# XXX_weight : [d1,d2,w1,w2] means:
# if XXX < d1, weight = w1
# if d1 <= XXX < d2, weight = w1 + (XXX-d1)/(d2-d1)*(w2-w1), (linear interpolation)
# if d2 <= XXX , weight = w2
# You can easily set w1 = w2 = 1.0 to normalize the weight related to XXX.
# -------------- using absolute traveltime data --------------
abs_time:
use_abs_time: true # 'true' for using absolute traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
residual_weight: [1, 3, 1, 1] # XXX is the absolute traveltime residual (second) = abs(t^{obs}_{n,i} - t^{syn}_{n,j})
distance_weight: [50, 150, 1, 1] # XXX is epicenter distance (km) between the source and receiver related to the data
# -------------- using common source differential traveltime data --------------
cs_dif_time:
use_cs_time: false # 'true' for using common source differential traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
residual_weight: [1, 3, 1, 0.1] # XXX is the common source differential traveltime residual (second) = abs(t^{obs}_{n,i} - t^{obs}_{n,j} - t^{syn}_{n,i} + t^{syn}_{n,j}).
azimuthal_weight: [15, 30, 1, 0.1] # XXX is the azimuth difference between two separate stations related to the common source.
# -------------- using common receiver differential traveltime data --------------
cr_dif_time:
use_cr_time: false # 'true' for using common receiver differential traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
residual_weight: [1, 3, 1, 0.1] # XXX is the common receiver differential traveltime residual (second) = abs(t^{obs}_{n,i} - t^{obs}_{m,i} - t^{syn}_{n,i} + t^{syn}_{m,i})
azimuthal_weight: [15, 30, 1, 0.1] # XXX is the azimuth difference between two separate sources related to the common receiver.
# -------------- global weight of different types of data (to balance the weight of different data) --------------
global_weight:
balance_data_weight: true # true: over the total weight of the each type of the data. false: use original weight (below weight for each type of data needs to be set)
abs_time_weight: 1 # weight of absolute traveltime data after balance, default: 1.0
cs_dif_time_local_weight: 1 # weight of common source differential traveltime data after balance, default: 1.0
cr_dif_time_local_weight: 1 # weight of common receiver differential traveltime data after balance, default: 1.0
teleseismic_weight: 1 # weight of teleseismic data after balance, default: 1.0 (exclude in this version)The parameters which will be optimized can be specified in this section.
# -------------- inversion parameters --------------
update_slowness : true # update slowness (velocity) or not. default: true
update_azi_ani : true # update azimuthal anisotropy (xi, eta) or not. default: false
# -------------- for teleseismic inversion (under development) --------------
# depth_taper : [d1,d2] means:
# if XXX < d1, kernel <- kernel * 0.0
# if d1 <= XXX < d2, kernel <- kernel * (XXX-d1)/(d2-d1), (linear interpolation)
# if d2 <= XXX , kernel <- kernel * 1.0
# You can easily set d1 = -200, d1 = -100 to remove this taper.
depth_taper : [-200, -100]Relocation parameters setting
The settings used for relocation is defined separately in this section.
relocation: # update earthquake hypocenter and origin time (when run_mode : 2 and 3)
min_Ndata: 4 # if the number of data of the earthquake is less than <min_Ndata>, the earthquake will not be relocated. defaut value: 4
# relocation_strategy
step_length : 0.01 # step length of relocation perturbation at each iteration. 0.01 means maximum 1% perturbation for each iteration.
step_length_decay : 0.9 # if objective function increase, step size -> step length * step_length_decay. default: 0.9
rescaling_dep_lat_lon_ortime : [10, 10, 10, 1] # The perturbation is related to <rescaling_dep_lat_lon_ortime>. Unit: km,km,km,second
max_change_dep_lat_lon_ortime : [5, 5, 5, 0.5] # the change of dep,lat,lon,ortime do not exceed max_change. Unit: km,km,km,second
max_iterations : 100 # maximum number of iterations for relocation
tol_gradient : 0.0001 # if the norm of gradient is smaller than the tolerance, the iteration of relocation terminates
# -------------- using absolute traveltime data --------------
abs_time:
use_abs_time : true # 'true' for using absolute traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
residual_weight : [1, 3, 1, 0.1] # XXX is the absolute traveltime residual (second) = abs(t^{obs}_{n,i} - t^{syn}_{n,j})
distance_weight : [50, 150, 1, 0.1] # XXX is epicenter distance (km) between the source and receiver related to the data
# -------------- using common receiver differential traveltime data --------------
cr_dif_time:
use_cr_time : true # 'true' for using common receiver differential traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
residual_weight : [1, 3, 1, 0.1] # XXX is the common receiver differential traveltime residual (second) = abs(t^{obs}_{n,i} - t^{obs}_{m,i} - t^{syn}_{n,i} + t^{syn}_{m,i})
azimuthal_weight : [10, 30, 1, 0.1] # XXX is the azimuth difference between two separate sources related to the common receiver.
# -------------- global weight of different types of data (to balance the weight of different data) --------------
global_weight:
balance_data_weight: false # true: over the total weight of the each type of the data. no: use original weight (below weight for each type of data needs to be set)
abs_time_local_weight: 1 # weight of absolute traveltime data for relocation after balance, default: 1.0
cr_dif_time_local_weight: 1 # weight of common receiver differential traveltime data for relocation after balance, default: 1.0Inversion strategy setting for tomography and relocation
TomoATT may also run the inversion and source position relocation iteratively. Those parameters are effective only with run_mode == 3.
inversion_strategy: # update model parameters and earthquake hypocenter iteratively (when run_mode : 3)
inv_mode : 0 # 0 for update model parameters and relocation iteratively. (other options for future work)
# for inv_mode : 0, parameters below are required
inv_mode_0: # update model for <model_update_N_iter> steps, then update location for <relocation_N_iter> steps, and repeat the process for <max_loop> loops.
model_update_N_iter : 1
relocation_N_iter : 1
max_loop : 10Eikonal solver setup
These are the parameters for the eikonal solver.
The first-order upwind scheme (stencil_type == 1) is only available with stencil_order == 1 and sweep_type == 0.
calculation:
convergence_tolerance: 0.0001 # threshold value for checking the convergence for each forward/adjoint run
max_iterations: 500 # number of maximum iteration for each forward/adjoint run
stencil_order: 3 # order of stencil, 1 or 3
stencil_type: 0 # 0: general sweeping scheme, 1: first-order upwind scheme (only sweep_type 0 is supported)
sweep_type: 1 # 0: legacy, 1: cuthill-mckee with shm parallelization2. source receiver file
Source receiver file is a file which defines source and receiver positions and arrival times.
Below is a brief example of the file format:
1 1992 1 1 2 43 56.900 1.8000 98.9000 137.00 2.80 8 305644
1 1 PCBI 1.8900 98.9253 1000.0000 P 18.000
1 2 MRPI 1.6125 99.3172 1100.0000 P 19.400
1 3 HUTI 2.3153 98.9711 1600.0000 P 19.200
....
The first line is the description of a seismic event. Then following lines are traveltime records for this event. The means of each column are:
Source line : id_src year month day hour min sec lat lon dep_km magnitude num_recs id_event (weight)
Receiver line : id_src id_rec name_rec lat lon elevation_m phase arrival_time_sec (weight)lon and lat are in degree.
num_recs (number of receivers for this source) need to be the same with the number of receiver lines.
The last column of both source and receiver line is for put weight (on objective function). These are the optional column and set to be 1.0 if not written.
name_rec need to be different for each station. Otherwise the later station will overwrite the previous description.
If the source position is out of the global domain (defined in input parameter file.), the code will flag this event as a teleseismic event and run the dedicated routine for teleseimic event. For teleseismic event, swap_src_rec will be ignored for this event (as the teleseismic case, a source is not a point but boundary surfaces).
For using double-difference traveltime, following format may be used:
0 1998 1 1 0 0 0.00 0.1000 0.5000 30.0000 3.00 3 s0
0 0 r0 0.1500 0.5000 300.0 P 7.161
0 1 r1 0.2900 0.5000 300.0 2 r2 0.4300 0.5000 300.0 P,cs -1.122
0 3 r3 0.5700 0.5000 300.0 1 s1 0.3667 0.5000 30.0000 P,cr 5.506The first line is the source definition.
The second line is the absolute traveltime without weight value.
The third line is the common source double-difference traveltime, which need to be defined with P,cs at the phase column.
The fourth line is the common receiver double-difference traveltime, which need to be defined with P,cr at the phase column.
3. initial model file
Initial model file is used for defining parameters of input mode.
Necessary parameters are vel (velocity in km/s), eta, xi.
inital model file in HDF5 format
In HDF5 I/O mode (output_file_format: 0 in input parameter file), all the necessary parameters should be saved in one single .h5 file, with exactly the same name of dataset with parameter name written above.
The dimension of dataset should be the same with n_rtp in input parameter file.
Please refer the examples/inversion_small/make_test_model.py for details.
initial model file in ASCII format
In ASCII I/O mode (output_file_format: 1 in input parameter file), all the necessary parameters should be save in one single ASCII file.
The number of rows in the file need to be equivalent with the number of global nodes (i.e. n_rtp[0]*n_rtp[1]*n_rtp[2]).
The node order should be:
# write nodes in rtp
for ir in range(n_rtp[0]): # number of nodes on r direction
for it in range(n_rtp[1]): # number of nodes on t direction
for ip in range(n_rtp[2]): # number of nodes on p direction
# write out parameters
# eta xi zeta fun fac_a fac_b fac_c fac_f
Complete example may be found examples/inversion_small_ASCII/make_test_model.py.
Output files
Calculated travel times at the stations will be writen in (source receiver file)_out.dat on the column for travel times.
Volumetric result data files are saved in OUTPUT_FILES directory.
As the node order in the output file is not in the global domain but each subdomains, it is necessary to do a small post-processing for extracting slices.
utils/tomoatt_data_retrieval.py includes functions for this post processing.
Please refer the concrete example in inversion_small/data_post_process.py for HDF5 mode, and inversion_small_ASCII/data_post_process.py for ASCII mode.
Only the final iteration result will be saved in 3D matrix form as final_model.h5 or final_model.dat, thus it is not necessary to do post-processing for the final result.
HDF5 I/O mode
In HDF5 mode, the code will carry out collective writing from all MPI processes into one single output file, which will try to maximize the I/O bandwidth for efficient I/O.
TomoATT produces output files like below:
- out_data_grid.h5 : grid coordinate and connectivity data
- out_data_sim.h5 : field data
- out_data_sim.xmf : XDMF index data for visualizing 3D data. This may be open by Paraview.
- final_model.h5 : final model parameters
Travel time field for i-th source may be visualized by reading OUTPUT_FILES/out_data_sim_i.xmf.
All the inversed parameters from all the sources and receivers are saved in out_data_sim_0.xmf.
Internal composition of .h5 data may be indicated by h5ls -r command.
The composition of out_data_grid.h5 is :
/ Group
/Mesh Group
/Mesh/elem_conn Dataset {21609, 9} # node connectivity used for visualization
/Mesh/node_coords_p Dataset {26010} # longitude [degree]
/Mesh/node_coords_r Dataset {26010} # radious [km]
/Mesh/node_coords_t Dataset {26010} # latiude [degree]
/Mesh/node_coords_x Dataset {26010} # xyz coordinate in WGS84
/Mesh/node_coords_y Dataset {26010}
/Mesh/node_coords_z Dataset {26010}
/Mesh/procid Dataset {26010} # mpi processor idout_data_sim.h5 includes the fields below :
/model Group
/model/eta_inv_000i Dataset {25000} # eta at i-th inversion
/model/xi_inv_000i Dataset {25000} # xi at i-th inversion
/model/vel_inv_000i Dataset {25000} # velocity field at i-th inversion
/model/Keta_inv_000i Dataset {25000} # kernel for eta at i-th inversion
/model/Kxi_inv_000i Dataset {25000} # kernel for xi at i-th inversion
/model/Ks_inv_000i Dataset {25000} # kernel for slowness field at i-th inversion
/model/Kxi_update_inv_000i Dataset {25000} # Smoothed kernel for xi at i-th inversion
/model/Keta_update_inv_000i Dataset {25000} # Smoothed kernel for eta at i-th inversion
/model/Ks_update_inv_000i Dataset {25000} # Smoothed kernel for slowness field at i-th inversion
/src_rec_j Group # data for j-th source
/src_rec_j/T_res_inv_000i Dataset {25000} # calculated traveltime field at i-th inversion
/src_rec_j/adjoint_field_inv_000i Dataset {25000} # adjoint field at i-th inversionthen final_model.h5 is :
eta Dataset {10, 50, 50}
vel Dataset {10, 50, 50}
xi Dataset {10, 50, 50}The internal composition of the final_model.h5 is exactly the same with the initial model file. Thus it is possible to use the final model file as the initial model file for the next run by specifying initial_model_file in input parameter file.
ASCII I/O mode
In ASCII mode, code will be do independent writing, (i.e. each MPI process do I/O process sequencially) in a single output file.
The files that TomoATT creates in OUPUT_FILES directory are :
out_grid_conn.dat # node connectivity
out_grid_ptr.dat # grid coordinate lon (degree), lat (degree), radius (km)
out_grid_xyz.dat # grid coordinate in WGS84
eta_inv_000i.dat # eta
xi_inv_000i.dat # xi
vel_inv_000i.dat # velocity