seisflows

plot code

fig1. Fig2.

%cd /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/specfem2d_workdir/
init = Model('OUTPUT_FILES_INIT')
init.plot2d('vs')
true = Model('OUTPUT_FILES_TRUE')
true.plot2d('vs')
#true = Model('OUTPUT_FILES')
#true.plot2d('vs')

fig3. fig4

%cd /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/output_1/
m_init = Model('MODEL_INIT')
m_true = Model('MODEL_TRUE')
#m_init.save('/home/zhangzhiyu/MyProjects/Seisflows/work/day_3/3/output/m_init.npz')
#m_true.plot2d('vs')

g_01 = Model('GRADIENT_01')
g_01.coordinates = {}
g_01.coordinates["x"] = m_init.coordinates["x"]
g_01.coordinates["z"] = m_init.coordinates["z"]
g_01.plot2d('vs_kernel')

m_01 = Model('MODEL_01')
m_01.coordinates = {}
m_01.coordinates["x"] = m_init.coordinates["x"]
m_01.coordinates["z"] = m_init.coordinates["z"]
m_01.plot2d('vs')

#m_01.save('/home/zhangzhiyu/MyProjects/Seisflows/work/day_3/3/output/m.npz')

parameters.yaml

# //////////////////////////////////////////////////////////////////////////////
#
#                        SeisFlows YAML Parameter File
#
# //////////////////////////////////////////////////////////////////////////////
#
# Modules correspond to the structure of the source code, and determine
# SeisFlows' behavior at runtime. Each module requires its own sub-parameters.
#
# .. rubric::
#   - Determine available options for modules by running:
#       > seisflows print modules
#   - Auto-fill with docstrings and default values (recommended) by running:
#       > seisflows configure
#   - Swap out module parameters for a configured parameter file by running:
#       > seisflows swap {module} {name} (e.g., seisflows swap solver specfem3d)
#   - To set values as NoneType, use: null
#   - To set values as infinity, use: inf
#
#                                    MODULES
#                                    ///////
# workflow (str):    The types and order of functions for running SeisFlows
# system (str):      Computer architecture of the system being used
# solver (str):      External numerical solver to use for waveform simulations
# preprocess (str):  Preprocessing schema for waveform data
# optimize (str):    Optimization algorithm for the inverse problem
# ==============================================================================
workflow: inversion
system: workstation
solver: specfem2d
preprocess: pyaflowa
optimize: LBFGS
# =============================================================================
#
#    Forward Workflow
#    ----------------
#    Defines foundational structure for Workflow module. When used standalone 
#    is in charge of running forward solver in parallel and (optionally) 
#    calculating data-synthetic misfit and adjoint sources.
#
#    Parameters
#    ----------
#    :type modules: list of module
#    :param modules: instantiated SeisFlows modules which should have been
#        generated by the function `seisflows.config.import_seisflows` with a
#        parameter file generated by seisflows.configure
#    :type data_case: str
#    :param data_case: How to address 'data' in the workflow, available options:
#        'data': real data will be provided by the user in
#        `path_data/{source_name}` in the same format that the solver will
#        produce synthetics (controlled by `solver.format`) OR
#        synthetic': 'data' will be generated as synthetic seismograms using
#        a target model provided in `path_model_true`. If None, workflow will
#        not attempt to generate data.
#    :type stop_after: str
#    :param stop_after: optional name of task in task list (use
#        `seisflows print tasks` to get task list for given workflow) to stop
#        workflow after, allowing user to prematurely stop a workflow to explore
#        intermediate results or debug.
#    :type export_traces: bool
#    :param export_traces: export all waveforms that are generated by the
#        external solver to `path_output`. If False, solver traces stored in
#        scratch may be discarded at any time in the workflow
#    :type export_residuals: bool
#    :param export_residuals: export all residuals (data-synthetic misfit) that
#        are generated by the external solver to `path_output`. If False,
#        residuals stored in scratch may be discarded at any time in the 
#        workflow
#
#        
#    Migration Workflow
#    ------------------
#    Run forward and adjoint solver to produce event-dependent misfit kernels.
#    Sum and postprocess kernels to produce gradient. In seismic exploration
#    this is 'reverse time migration'.
#
#    Parameters
#    ----------
#    :type export_gradient: bool
#    :param export_gradient: export the gradient after it has been generated
#        in the scratch directory. If False, gradient can be discarded from
#        scratch at any time in the workflow
#    :type export_kernels: bool
#    :param export_kernels: export each sources event kernels after they have
#        been generated in the scratch directory. If False, gradient can be
#        discarded from scratch at any time in the workflow
#
#        
#    Inversion Workflow
#    ------------------
#    Peforms iterative nonlinear inversion using the machinery of the Forward
#    and Migration workflows, as well as a built-in optimization library.
#
#    Parameters
#    ----------
#    :type start: int
#    :param start: start inversion workflow at this iteration. 1 <= start <= inf
#    :type end: int
#    :param end: end inversion workflow at this iteration. start <= end <= inf
#    :type iteration: int
#    :param iteration: The current iteration of the workflow. If NoneType, takes
#        the value of `start` (i.e., first iteration of the workflow). User can
#        also set between `start` and `end` to resume a failed workflow.
#    :type thrifty: bool
#    :param thrifty: a thrifty inversion skips the costly intialization step
#        (i.e., forward simulations and misfit quantification) if the final
#        forward simulations from the previous iterations line search can be
#        used in the current one. Requires L-BFGS optimization.
#    :type export_model: bool
#    :param export_model: export best-fitting model from the line search to disk.
#        If False, new models can be discarded from scratch at any time.
#
#        
# =============================================================================
data_case: synthetic
stop_after: null
export_traces: False
export_residuals: False
export_gradient: True
export_kernels: False
start: 1
end: 30
export_model: True
thrifty: False
iteration: 1
# =============================================================================
#
#    Workstation System
#    ------------------
#    Defines foundational structure for System module. When used standalone,
#    runs solver tasks either in serial (if `nproc`==1; i.e., without MPI) or in
#    parallel (if `nproc`>1; i.e., with MPI). All other tasks are run in serial.
#
#    Parameters
#    ----------
#    :type ntask: int
#    :param ntask: number of individual tasks/events to run during workflow.
#        Must be <= the number of source files in `path_specfem_data`
#    :type nproc: int
#    :param nproc: number of processors to use for each simulation. Choose 1 for
#        serial simulations, and `nproc`>1 for parallel simulations.
#    :type mpiexec: str
#    :param mpiexec: MPI executable on system. Defaults to 'mpirun -n ${NPROC}'
#    :type log_level: str
#    :param log_level: logger level to pass to logging module.
#        Available: 'debug', 'info', 'warning', 'critical'
#    :type verbose: bool
#    :param verbose: if True, formats the log messages to include the file
#        name, line number and message type. Useful for debugging but
#        also very verbose.
#
#        
# =============================================================================
ntask: 4
nproc: 4
mpiexec: mpirun
log_level: DEBUG
verbose: False
# =============================================================================
#
#    Solver SPECFEM
#    --------------
#    Defines foundational structure for Specfem-based solver module. 
#    Generalized SPECFEM interface to manipulate SPECFEM2D/3D/3D_GLOBE w/ Python
#
#    Parameters
#    ----------
#    :type data_format: str
#    :param data_format: data format for reading traces into memory.
#        Available: ['SU': seismic unix format, 'ASCII': human-readable ascii]
#    :type materials: str
#    :param materials: Material parameters used to define model. Available:
#        ['ELASTIC': Vp, Vs, 'ACOUSTIC': Vp, 'ISOTROPIC', 'ANISOTROPIC']
#    :type density: bool
#    :param density: How to treat density during inversion. If True, updates
#        density during inversion. If False, keeps it constant.
#        TODO allow density scaling during an inversion
#    :type attenuation: bool
#    :param attenuation: How to treat attenuation during inversion.
#        if True, turns on attenuation during forward simulations only. If
#        False, attenuation is always set to False. Requires underlying
#        attenution (Q_mu, Q_kappa) model
#    :type smooth_h: float
#    :param smooth_h: Gaussian half-width for horizontal smoothing in units
#        of meters. If 0., no smoothing applied. Only applicable for workflows:
#        ['migration', 'inversion'], ignored for 'forward' workflow.
#    :type smooth_h: float
#    :param smooth_v: Gaussian half-width for vertical smoothing in units
#        of meters. Only applicable for workflows: ['migration', 'inversion'],
#        ignored for 'forward' workflow.
#    :type components: str
#    :param components: components to consider and tag data with. Should be
#        string of letters such as 'RTZ'
#    :type solver_io: str
#    :param solver_io: format of model/kernel/gradient files expected by the
#        numerical solver. Available: ['fortran_binary': default .bin files].
#        TODO: ['adios': ADIOS formatted files]
#    :type source_prefix: str
#    :param source_prefix: prefix of source/event/earthquake files. If None,
#        will attempt to guess based on the specific solver chosen.
#    :type mpiexec: str
#    :param mpiexec: MPI executable used to run parallel processes. Should also
#        be defined for the system module
#
#        
#    Solver SPECFEM2D
#    ----------------
#    SPECFEM2D-specific alterations to the base SPECFEM module
#
#    Parameters
#    ----------
#    :type source_prefix: str
#    :param source_prefix: Prefix of source files in path SPECFEM_DATA. Defaults
#        to 'SOURCE'
#    :type multiples: bool
#    :param multiples: set an absorbing top-boundary condition
#
#        
# =============================================================================
data_format: ascii
materials: elastic
density: False
attenuation: False
smooth_h: 40000.0
smooth_v: 40000.0
components: XZ
source_prefix: SOURCE
multiples: False
# =============================================================================
#
#    Pyaflowa Preprocess
#    -------------------
#    Preprocessing and misfit quantification using Python's Adjoint Tomography
#    Operations Assistance (Pyatoa)
#
#    Parameters
#    ----------
#    :type min_period: float
#    :param min_period: Minimum filter corner in unit seconds. Bandpass
#    filter if set with `max_period`, highpass filter if set without
#    `max_period`, no filtering if not set and `max_period also not set
#    :type max_period: float
#    :param max_period: Maximum filter corner in unit seconds. Bandpass
#        filter if set with `min_period`, lowpass filter if set without
#        `min_period`, no filtering if not set and `min_period also not set
#   :type filter_corners: int
#    :param filter_corners: number of filter corners applied to filtering
#    :type client: str
#    :param client: Client name for ObsPy FDSN data gathering. Pyatoa will
#        attempt to collect waveform and metadata based on network and
#        station codes provided in the SPECFEM STATIONS file. If set None,
#        no FDSN gathering will be attempted
#    :type rotate: bool
#    :param rotate: Attempt to rotate waveform components from NEZ -> RTZ
#    :type pyflex_preset: str
#    :param pyflex_preset: Parameter map for misfit window configuration
#        defined by Pyflex. IF None, misfit and adjoint sources will be
#        calculated on whole traces. For available choices, see Pyatoa docs
#        page (pyatoa.rtfd.io)
#    :type fix_windows: bool or str
#    :param fix_windows: How to address misfit window evaluation at each
#        evaluation. Options to re-use misfit windows collected during an
#        inversion, available options:
#        [True, False, 'ITER', 'ONCE']
#        True: Re-use windows after first evaluation (i01s00);
#        False: Calculate new windows each evaluation;
#        'ITER': Calculate new windows at first evaluation of
#        each iteration (e.g., i01s00... i02s00...
#        'ONCE': Calculate new windows at first evaluation of
#        the workflow, i.e., at self.par.BEGIN
#    :type adj_src_type: str
#    :param adj_src_type: Adjoint source type to evaluate misfit, defined by
#        Pyadjoint. Currently available options: ['cc': cross-correlation,
#        'mt': multitaper, 'wav': waveform']
#    :type plot: bool
#    :param plot: plot waveform figures and source receiver maps during
#        the preprocessing stage
#    :type pyatoa_log_level: str
#    :param pyatoa_log_level: Log level to set Pyatoa, Pyflex, Pyadjoint.
#        Available: ['null': no logging, 'warning': warnings only,
#        'info': task tracking, 'debug': log all small details (recommended)]
#    :type unit_output: str
#    :param unit_output: Data units. Must match the synthetic output of
#        external solver. Available: ['DISP': displacement, 'VEL': velocity,
#        'ACC': acceleration]
#    :type export_datasets: bool
#    :param export_datasets: periodically save the output ASDFDataSets which
#        contain data, metadata and results collected during the
#        preprocessing procedure
#    :type export_figures: bool
#    :param export_figures: periodically save the output basemaps and
#        data-synthetic waveform comparison figures
#    :type export_log_files: bool
#    :param export_log_files: periodically save log files created by Pyatoa
#
#        
# =============================================================================
min_period: 10.0
max_period: 200.0
filter_corners: 4
client: null
rotate: False
pyflex_preset: null
fix_windows: False
adj_src_type: waveform
plot: True
pyatoa_log_level: DEBUG
unit_output: DISP
max_workers: null
export_datasets: True
export_figures: True
export_log_files: True
# =============================================================================
#
#    Gradient Optimization
#    ---------------------
#    Defines foundational structure for Optimization module. Applies a 
#    gradient/steepest descent optimization algorithm.
#
#    Parameters
#    ----------
#    :type line_search_method: str
#    :param line_search_method: chosen line_search algorithm. Currently available
#        are 'bracket' and 'backtrack'. See seisflows.plugins.line_search
#        for all available options
#    :type preconditioner: str
#    :param preconditioner: algorithm for preconditioning gradients. Currently
#        available: 'diagonal'. Requires `path_preconditioner` to point to a
#        set of files that define the preconditioner, formatted the same as the
#        input model
#    :type step_count_max: int
#    :param step_count_max: maximum number of trial steps to perform during
#        the line search before a change in line search behavior is
#        considered, or a line search is considered to have failed.
#    :type step_len_init: float
#    :param step_len_init: initial line search step length guess, provided
#        as a fraction of current model parameters.
#    :type step_len_max: float
#    :param step_len_max: maximum allowable step length during the line
#        search. Set as a fraction of the current model parameters
#
#        
#    L-BFGS Optimization
#    -------------------
#    Limited memory BFGS nonlienar optimization algorithm
#
#    Parameters
#    ----------
#    :type lbfgs_mem: int
#    :param lbfgs_mem: L-BFGS memory. Max number of previous gradients to
#        retain in local memory for approximating the objective function.
#    :type lbfgs_max: L-BFGS periodic restart interval. Must be
#        1 <= lbfgs_max <= infinity.
#    :type lbfgs_thresh: L-BFGS angle restart threshold. If the angle between
#        the current and previous search direction exceeds this value,
#        optimization algorithm will be restarted.
#
#        
# =============================================================================
preconditioner: null
step_count_max: 5
step_len_init: 0.05
step_len_max: 0.5
line_search_method: Backtrack
LBFGS_mem: 3
LBFGS_max: inf
LBFGS_thresh: 0.0
# =============================================================================
#
#	 Paths
#	 -----
#    :type workdir: str
#    :param workdir: working directory in which to perform a SeisFlows workflow.
#        SeisFlows internal directory structure will be created here. Default cwd
#    :type path_output: str
#    :param path_output: path to directory used for permanent storage on disk.
#        Results and exported scratch files are saved here.
#    :type path_data: str
#    :param path_data: path to any externally stored data required by the solver
#    :type path_state_file: str
#    :param path_state_file: path to a text file used to track the current
#        status of a workflow (i.e., what functions have already been completed),
#        used for checkpointing and resuming workflows
#    :type path_model_init: str
#    :param path_model_init: path to the starting model used to calculate the
#        initial misfit. Must match the expected `solver_io` format.
#    :type path_model_true: str
#    :param path_model_true: path to a target model if `case`=='synthetic' and
#        a set of synthetic 'observations' are required for workflow.
#    :type path_eval_grad: str
#    :param path_eval_grad: scratch path to store files for gradient evaluation,
#        including models, kernels, gradient and residuals.
#        :type path_mask: str
#    :param path_mask: optional path to a masking function which is used to
#        mask out or scale parts of the gradient. The user-defined mask must
#        match the file format of the input model (e.g., .bin files).
#        :type path_eval_func: str
#    :param path_eval_func: scratch path to store files for line search objective
#        function evaluations, including models, misfit and residuals
#        
#    :type path_output_log: str
#    :param path_output_log: path to a text file used to store the outputs of
#        the package wide logger, which are also written to stdout
#    :type path_par_file: str
#    :param path_par_file: path to parameter file which is used to instantiate
#        the package
#    :type path_log_files: str
#    :param path_log_files: path to a directory where individual log files are
#        saved whenever a number of parallel tasks are run on the system.
#        
#    :type path_data: str
#    :param path_data: path to any externally stored data required by the solver
#    :type path_specfem_bin: str
#    :param path_specfem_bin: path to SPECFEM bin/ directory which
#        contains binary executables for running SPECFEM
#    :type path_specfem_data: str
#    :param path_specfem_data: path to SPECFEM DATA/ directory which must
#        contain the CMTSOLUTION, STATIONS and Par_file files used for
#        running SPECFEM
#            
#    :type path_preprocess: str
#    :param path_preprocess: scratch path for preprocessing related steps
#        
#    :type path_preconditioner: str
#    :param path_preconditioner: optional path to a set of preconditioner files
#        formatted the same as the input model (or output model of solver).
#        Required to exist and contain files if `preconditioner`==True
#            
# =============================================================================
path_workdir: /home/erbiaoger/MyProjects/Seisflows/work/day_3/4
path_scratch: /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/scratch
path_eval_grad: /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/scratch/eval_grad
path_output: /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/output
path_model_init: /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/specfem2d_workdir/OUTPUT_FILES_INIT
path_model_true: /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/specfem2d_workdir/OUTPUT_FILES_TRUE
path_state_file: /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/sfstate.txt
path_data: null
path_mask: null
path_eval_func: /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/scratch/eval_func
path_par_file: /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/parameters.yaml
path_log_files: /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/logs
path_output_log: /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/sflog.txt
path_specfem_bin: /home/erbiaoger/MyProjects/Seisflows/specfem2d/bin
path_specfem_data: /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/specfem2d_workdir/DATA
path_solver: /home/erbiaoger/MyProjects/Seisflows/work/day_3/44/scratch/solver
path_preconditioner: null

Par_file

#-----------------------------------------------------------
#
# Simulation input parameters
#
#-----------------------------------------------------------

# title of job
title                           = Tape-Liu-Tromp (GJI 2007)

# forward or adjoint simulation
# 1 = forward, 2 = adjoint, 3 = both simultaneously
# note: 2 is purposely UNUSED (for compatibility with the numbering of our 3D codes)
SIMULATION_TYPE                 = 1
# 0 = regular wave propagation simulation, 1/2/3 = noise simulation
NOISE_TOMOGRAPHY                = 0
# save the last frame, needed for adjoint simulation
SAVE_FORWARD                    = .false.

# parameters concerning partitioning
NPROC                           = 4              # number of processes

# time step parameters
# total number of time steps
NSTEP                           = 14000

# duration of a time step (see section "How to choose the time step" of the manual for how to do this)
DT                              = 2.5d-2

# time stepping
# 1 = Newmark (2nd order), 2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta), 3 = classical RK4 4th-order 4-stage Runge-Kutta
time_stepping_scheme            = 1

# set the type of calculation (P-SV or SH/membrane waves)
P_SV                            = .true.

# axisymmetric (2.5D) or Cartesian planar (2D) simulation
AXISYM                          = .false.

#-----------------------------------------------------------
#
# Mesh
#
#-----------------------------------------------------------

# Partitioning algorithm for decompose_mesh
PARTITIONING_TYPE               = 3              # SCOTCH = 3, ascending order (very bad idea) = 1

# number of control nodes per element (4 or 9)
NGNOD                           = 4

# creates/reads a binary database that allows to skip all time consuming setup steps in initialization
# 0 = does not read/create database
# 1 = creates database
# 2 = reads database
setup_with_binary_database      = 1

# available models
#   default       - define model using nbmodels below
#   ascii         - read model from ascii database file
#   binary        - read model from binary databse file
#   binary_voigt  - read Voigt model from binary database file
#   external      - define model using define_external_model subroutine
#   gll           - read GLL model from binary database file
#   legacy        - read model from model_velocity.dat_input
MODEL                           = gll

# Output the model with the requested type, does not save if turn to default or .false.
# (available output formats: ascii,binary,gll,legacy)
SAVE_MODEL                      = binary


#-----------------------------------------------------------
#
# Attenuation
#
#-----------------------------------------------------------

# attenuation parameters
ATTENUATION_VISCOELASTIC        = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
ATTENUATION_VISCOACOUSTIC       = .false.        # turn attenuation (viscoacousticity) on or off for non-poroelastic fluid parts of the model

# for viscoelastic or viscoacoustic attenuation
N_SLS                           = 3              # number of standard linear solids for attenuation (3 is usually the minimum)
ATTENUATION_f0_REFERENCE        = 5.196152422706633 # in case of attenuation, reference frequency in Hz at which the velocity values in the velocity model are given (unused otherwise); relevant only if source is a Dirac or a Heaviside, otherwise it is automatically set to f0 the dominant frequency of the source in the DATA/SOURCE file
READ_VELOCITIES_AT_f0           = .false.        # read seismic velocities at ATTENUATION_f0_REFERENCE instead of at infinite frequency (see user manual for more information)
USE_SOLVOPT                     = .false.        # use more precise but much more expensive way of determining the Q factor relaxation times, as in https://doi.org/10.1093/gji/ggw024

# for poroelastic attenuation
ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
Q0_poroelastic                  = 1              # quality factor for viscous attenuation (ignore it if you are not using a poroelastic material)
freq0_poroelastic               = 10             # frequency for viscous attenuation (ignore it if you are not using a poroelastic material)

# to undo attenuation and/or PMLs for sensitivity kernel calculations or forward runs with SAVE_FORWARD
# use the flag below. It performs undoing of attenuation and/or of PMLs in an exact way for sensitivity kernel calculations
# but requires disk space for temporary storage, and uses a significant amount of memory used as buffers for temporary storage.
# When that option is on the second parameter indicates how often the code dumps restart files to disk (if in doubt, use something between 100 and 1000).
UNDO_ATTENUATION_AND_OR_PML     = .false.
NT_DUMP_ATTENUATION             = 500

# Instead of reconstructing the forward wavefield, this option reads it from the disk using asynchronous I/O.
# Outperforms conventional mode using a value of NTSTEP_BETWEEN_COMPUTE_KERNELS high enough.
NO_BACKWARD_RECONSTRUCTION      = .false.

#-----------------------------------------------------------
#
# Sources
#
#-----------------------------------------------------------

# source parameters
NSOURCES                        = 1              # number of sources (source information is then read from the DATA/SOURCE file)
force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)

# use an existing initial wave field as source or start from zero (medium initially at rest)
initialfield                    = .false.
add_Bielak_conditions_bottom    = .false.        # add Bielak conditions or not if initial plane wave
add_Bielak_conditions_right     = .false.
add_Bielak_conditions_top       = .false.
add_Bielak_conditions_left      = .false.

# acoustic forcing
ACOUSTIC_FORCING                = .false.        # acoustic forcing of an acoustic medium with a rigid interface

# noise simulations - type of noise source time function:
# 0=external (S_squared), 1=Ricker(second derivative), 2=Ricker(first derivative), 3=Gaussian, 4=Figure 2a of Tromp et al. 2010
# (default value 4 is chosen to reproduce the time function from Fig 2a of "Tromp et al., 2010, Noise Cross-Correlation Sensitivity Kernels")
noise_source_time_function_type = 4

# moving sources
# Set write_moving_sources_database to .true. if the generation of moving source databases takes
# a long time. Then the simulation is done in two steps: first you run the code and it writes the databases to file
# (in DATA folder by default). Then you rerun the code and it will read the databases in there directly possibly
# saving a lot of time.
# This is only useful for GPU version (for now)
write_moving_sources_database   = .false.

#-----------------------------------------------------------
#
# Receivers
#
#-----------------------------------------------------------

# receiver set parameters for recording stations (i.e. recording points)
# seismotype : record 1=displ 2=veloc 3=accel 4=pressure 5=curl of displ 6=the fluid potential
seismotype                      = 1              # several values can be chosen. For example : 1,2,4

# interval in time steps for writing of seismograms
# every how many time steps we save the seismograms
# (costly, do not use a very small value; if you use a very large value that is larger than the total number
#  of time steps of the run, the seismograms will automatically be saved once at the end of the run anyway)
NTSTEP_BETWEEN_OUTPUT_SEISMOS   = 100

# set to n to reduce the sampling rate of output seismograms by a factor of n
# defaults to 1, which means no down-sampling
NTSTEP_BETWEEN_OUTPUT_SAMPLE    = 1

# so far, this option can only be used if all the receivers are in acoustic elements
USE_TRICK_FOR_BETTER_PRESSURE   = .false.

# use this t0 as earliest starting time rather than the automatically calculated one
USER_T0                         = 90.0d0

# seismogram formats
save_ASCII_seismograms          = .true.         # save seismograms in ASCII format or not
save_binary_seismograms_single  = .false.        # save seismograms in single precision binary format or not (can be used jointly with ASCII above to save both)
save_binary_seismograms_double  = .false.        # save seismograms in double precision binary format or not (can be used jointly with both flags above to save all)
SU_FORMAT                       = .false.        # output single precision binary seismograms in Seismic Unix format (adjoint traces will be read in the same format)

# use an existing STATION file found in ./DATA or create a new one from the receiver positions below in this Par_file
use_existing_STATIONS           = .false.

# number of receiver sets (i.e. number of receiver lines to create below)
nreceiversets                   = 1

# orientation
anglerec                        = 0.d0           # angle to rotate components at receivers
rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)

# first receiver set (repeat these 6 lines and adjust nreceiversets accordingly)
nrec                            = 40              # number of receivers
xdeb                            = 10000.      # first receiver x in meters
zdeb                            = 180000.      # first receiver z in meters
xfin                            = 470000.        # last receiver x in meters (ignored if only one receiver)
zfin                            = 180000.         # last receiver z in meters (ignored if only one receiver)
record_at_surface_same_vertical = .true.        # receivers inside the medium or at the surface


#-----------------------------------------------------------
#
# adjoint kernel outputs
#
#-----------------------------------------------------------

# save sensitivity kernels in ASCII format (much bigger files, but compatible with current GMT scripts) or in binary format
save_ASCII_kernels              = .false.

# since the accuracy of kernel integration may not need to respect the CFL, this option permits to save computing time, and memory with UNDO_ATTENUATION_AND_OR_PML mode
NTSTEP_BETWEEN_COMPUTE_KERNELS  = 1

# outputs approximate Hessian for preconditioning
APPROXIMATE_HESS_KL             = .false.

#-----------------------------------------------------------
#
# Boundary conditions
#
#-----------------------------------------------------------

# Perfectly Matched Layer (PML) boundaries
# absorbing boundary active or not
PML_BOUNDARY_CONDITIONS         = .false.
NELEM_PML_THICKNESS             = 3
ROTATE_PML_ACTIVATE             = .false.
ROTATE_PML_ANGLE                = 30.
# change the four parameters below only if you know what you are doing; they change the damping profiles inside the PMLs
K_MIN_PML                       = 1.0d0          # from Gedney page 8.11
K_MAX_PML                       = 1.0d0
damping_change_factor_acoustic  = 0.5d0
damping_change_factor_elastic   = 1.0d0
# set the parameter below to .false. unless you know what you are doing; this implements automatic adjustment of the PML parameters for elongated models.
# The goal is to improve the absorbing efficiency of PML for waves with large incidence angles, but this can lead to artefacts.
# In particular, this option is efficient only when the number of sources NSOURCES is equal to one.
PML_PARAMETER_ADJUSTMENT        = .false.

# Stacey ABC
STACEY_ABSORBING_CONDITIONS     = .true.

# periodic boundaries
ADD_PERIODIC_CONDITIONS         = .false.
PERIODIC_HORIZ_DIST             = 0.3597d0

#-----------------------------------------------------------
#
# Velocity and density models
#
#-----------------------------------------------------------

# number of model materials
nbmodels                        = 2
# available material types (see user manual for more information)
#   acoustic:              model_number 1 rho Vp 0  0 0 QKappa 9999 0 0 0 0 0 0 (for QKappa use 9999 to ignore it)
#   elastic:               model_number 1 rho Vp Vs 0 0 QKappa Qmu  0 0 0 0 0 0 (for QKappa and Qmu use 9999 to ignore them)
#   anisotropic:           model_number 2 rho c11 c13 c15 c33 c35 c55 c12 c23 c25   0 QKappa Qmu
#   anisotropic in AXISYM: model_number 2 rho c11 c13 c15 c33 c35 c55 c12 c23 c25 c22 QKappa Qmu
#   poroelastic:           model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu
#   tomo:                  model_number -1 0 0 A 0 0 0 0 0 0 0 0 0 0
#
# note: When viscoelasticity or viscoacousticity is turned on,
#       the Vp and Vs values that are read here are the UNRELAXED ones i.e. the values at infinite frequency
#       unless the READ_VELOCITIES_AT_f0 parameter above is set to true, in which case they are the values at frequency f0.
#
#       Please also note that Qmu is always equal to Qs, but Qkappa is in general not equal to Qp.
#       To convert one to the other see doc/Qkappa_Qmu_versus_Qp_Qs_relationship_in_2D_plane_strain.pdf and
#       utils/attenuation/conversion_from_Qkappa_Qmu_to_Qp_Qs_from_Dahlen_Tromp_959_960.f90.
1 1 2600.d0 5800.d0 3500.0d0 0 0 10.d0 10.d0 0 0 0 0 0 0
2 1 2600.d0 5800.d0 3500.0d0 0 0 10.d0 10.d0 0 0 0 0 0 0
#1 1 2600.d0 5600.d0 3200.0d0 0 0 10.d0 10.d0 0 0 0 0 0 0
#2 1 2600.d0 6000.d0 3800.0d0 0 0 10.d0 10.d0 0 0 0 0 0 0

# external tomography file
TOMOGRAPHY_FILE                 = ./DATA/tomo_file.xyz

# use an external mesh created by an external meshing tool or use the internal mesher
read_external_mesh              = .false.

#-----------------------------------------------------------
#
# PARAMETERS FOR EXTERNAL MESHING
##
#-----------------------------------------------------------

# data concerning mesh, when generated using third-party app (more info in README)
# (see also absorbing_conditions above)
mesh_file                       = ./DATA/ice_water_rock_1D/ice_water_rock_1D_elements   # file containing the mesh
nodes_coords_file               = ./DATA/ice_water_rock_1D/ice_water_rock_1D_nodes   # file containing the nodes coordinates
materials_file                  = ./DATA/ice_water_rock_1D/ice_water_rock_1D_material   # file containing the material number for each element
free_surface_file               = ./DATA/ice_water_rock_1D/ice_water_rock_1D_surface_free   # file containing the free surface
axial_elements_file             = ./DATA/axial_elements_file   # file containing the axial elements if AXISYM is true
absorbing_surface_file          = ./DATA/ice_water_rock_1D/ice_water_rock_1D_surface_absorb   # file containing the absorbing surface
acoustic_forcing_surface_file   = ./DATA/MSH/Surf_acforcing_Bottom_enforcing_mesh   # file containing the acoustic forcing surface
absorbing_cpml_file             = ./DATA/absorbing_cpml_file   # file containing the CPML element numbers
tangential_detection_curve_file = ./DATA/courbe_eros_nodes  # file containing the curve delimiting the velocity model

#-----------------------------------------------------------
#
# PARAMETERS FOR INTERNAL MESHING
#
#-----------------------------------------------------------

# file containing interfaces for internal mesh
interfacesfile                  = interfaces_Tape2007.dat

# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin                            = 0.d0           # abscissa of left side of the model
xmax                            = 480000.d0      # abscissa of right side of the model
nx                              = 40             # number of elements along X

# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom                    = .false.
absorbright                     = .true.
absorbtop                       = .true.
absorbleft                      = .true.

# define the different regions of the model in the (nx,nz) spectral-element mesh
nbregions                       = 8              # then set below the different regions and model number for each region
# format of each line: nxmin nxmax nzmin nzmax material_number
1 10  1 20 1
11 20  1 20 2
21 30  1 20 1
31 40  1 20 2
1 10  21 40 2
11 20  21 40 1
21 30  21 40 2
31 40  21 40 1

#-----------------------------------------------------------
#
# Display parameters
#
#-----------------------------------------------------------

# interval at which we output time step info and max of norm of displacement
# (every how many time steps we display information about the simulation. costly, do not use a very small value)
NTSTEP_BETWEEN_OUTPUT_INFO      = 400

# meshing output
output_grid_Gnuplot             = .false.        # generate a GNUPLOT file containing the grid, and a script to plot it
output_grid_ASCII               = .false.        # dump the grid in an ASCII text file consisting of a set of X,Y,Z points or not

# to plot total energy curves, for instance to monitor how CPML absorbing layers behave;
# should be turned OFF in most cases because a bit expensive
OUTPUT_ENERGY                   = .false.

# every how many time steps we compute energy (which is a bit expensive to compute)
NTSTEP_BETWEEN_OUTPUT_ENERGY    = 10

# Compute the field int_0^t v^2 dt for a set of GLL points and write it to file. Use
# the script utils/visualisation/plotIntegratedEnergyFile.py to watch. It is refreshed at the same time than the seismograms
COMPUTE_INTEGRATED_ENERGY_FIELD = .false.

#-----------------------------------------------------------
#
# Movies/images/snaphots visualizations
#
#-----------------------------------------------------------

# every how many time steps we draw JPEG or PostScript pictures of the simulation
# and/or we dump results of the simulation as ASCII or binary files (costly, do not use a very small value)
NTSTEP_BETWEEN_OUTPUT_IMAGES    = 200

# minimum amplitude kept in % for the JPEG and PostScript snapshots; amplitudes below that are muted
cutsnaps                        = 1.

#### for JPEG color images ####
output_color_image              = .true.         # output JPEG color image of the results every NTSTEP_BETWEEN_OUTPUT_IMAGES time steps or not
imagetype_JPEG                  = 3              # display 1=displ_Ux 2=displ_Uz 3=displ_norm 4=veloc_Vx 5=veloc_Vz 6=veloc_norm 7=accel_Ax 8=accel_Az 9=accel_norm 10=pressure
factor_subsample_image          = 1.0d0          # (double precision) factor to subsample or oversample (if set to e.g. 0.5) color images output by the code (useful for very large models, or to get nicer looking denser pictures)
USE_CONSTANT_MAX_AMPLITUDE      = .false.        # by default the code normalizes each image independently to its maximum; use this option to use the global maximum below instead
CONSTANT_MAX_AMPLITUDE_TO_USE   = 1.17d4         # constant maximum amplitude to use for all color images if the above USE_CONSTANT_MAX_AMPLITUDE option is true
POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in JPEG color images
DRAW_SOURCES_AND_RECEIVERS      = .true.         # display sources as orange crosses and receivers as green squares in JPEG images or not
DRAW_WATER_IN_BLUE              = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water in the case of ocean acoustics or in the case of offshore oil industry experiments (if off, display them as greyscale, as for elastic or poroelastic elements, for instance for acoustic-only oil industry models of solid media)
USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step (for instance to create movies in an easier way later)

#### for PostScript snapshots ####
output_postscript_snapshot      = .false.        # output Postscript snapshot of the results every NTSTEP_BETWEEN_OUTPUT_IMAGES time steps or not
imagetype_postscript            = 1              # display 1=displ vector 2=veloc vector 3=accel vector; small arrows are displayed for the vectors
meshvect                        = .false.        # display mesh on PostScript plots or not
modelvect                       = .false.        # display velocity model on PostScript plots or not
boundvect                       = .true.         # display boundary conditions on PostScript plots or not
interpol                        = .true.         # interpolation of the PostScript display on a regular grid inside each spectral element, or use the non-evenly spaced GLL points
pointsdisp                      = 6              # number of points in each direction for interpolation of PostScript snapshots (set to 1 for lower-left corner only)
subsamp_postscript              = 1              # subsampling of background velocity model in PostScript snapshots
sizemax_arrows                  = 1.d0           # maximum size of arrows on PostScript plots in centimeters
US_LETTER                       = .false.        # use US letter or European A4 paper for PostScript plots

#### for wavefield dumps ####
output_wavefield_dumps          = .false.        # output wave field to a text file (creates very big files)
imagetype_wavefield_dumps       = 1              # display 1=displ vector 2=veloc vector 3=accel vector 4=pressure
use_binary_for_wavefield_dumps  = .false.        # use ASCII or single-precision binary format for the wave field dumps

#-----------------------------------------------------------

# Ability to run several calculations (several earthquakes)
# in an embarrassingly-parallel fashion from within the same run;
# this can be useful when using a very large supercomputer to compute
# many earthquakes in a catalog, in which case it can be better from
# a batch job submission point of view to start fewer and much larger jobs,
# each of them computing several earthquakes in parallel.
# To turn that option on, set parameter NUMBER_OF_SIMULTANEOUS_RUNS to a value greater than 1.
# To implement that, we create NUMBER_OF_SIMULTANEOUS_RUNS MPI sub-communicators,
# each of them being labeled "my_local_mpi_comm_world", and we use them
# in all the routines in "src/shared/parallel.f90", except in MPI_ABORT() because in that case
# we need to kill the entire run.
# When that option is on, of course the number of processor cores used to start
# the code in the batch system must be a multiple of NUMBER_OF_SIMULTANEOUS_RUNS,
# all the individual runs must use the same number of processor cores,
# which as usual is NPROC in the Par_file,
# and thus the total number of processor cores to request from the batch system
# should be NUMBER_OF_SIMULTANEOUS_RUNS * NPROC.
# All the runs to perform must be placed in directories called run0001, run0002, run0003 and so on
# (with exactly four digits).
#
# Imagine you have 10 independent calculations to do, each of them on 100 cores; you have three options:
#
# 1/ submit 10 jobs to the batch system
#
# 2/ submit a single job on 1000 cores to the batch, and in that script create a sub-array of jobs to start 10 jobs,
# each running on 100 cores (see e.g. http://www.schedmd.com/slurmdocs/job_array.html )
#
# 3/ submit a single job on 1000 cores to the batch, start SPECFEM2D on 1000 cores, create 10 sub-communicators,
# cd into one of 10 subdirectories (called e.g. run0001, run0002,... run0010) depending on the sub-communicator
# your MPI rank belongs to, and run normally on 100 cores using that sub-communicator.
#
# The option below implements 3/.
#
NUMBER_OF_SIMULTANEOUS_RUNS     = 1

# if we perform simultaneous runs in parallel, if only the source and receivers vary between these runs
# but not the mesh nor the model (velocity and density) then we can also read the mesh and model files
# from a single run in the beginning and broadcast them to all the others; for a large number of simultaneous
# runs for instance when solving inverse problems iteratively this can DRASTICALLY reduce I/Os to disk in the solver
# (by a factor equal to NUMBER_OF_SIMULTANEOUS_RUNS), and reducing I/Os is crucial in the case of huge runs.
# Thus, always set this option to .true. if the mesh and the model are the same for all simultaneous runs.
# In that case there is no need to duplicate the mesh and model file database (the content of the DATABASES_MPI
# directories) in each of the run0001, run0002,... directories, it is sufficient to have one in run0001
# and the code will broadcast it to the others)
BROADCAST_SAME_MESH_AND_MODEL   = .true.

#-----------------------------------------------------------

# set to true to use GPUs
GPU_MODE                        = .false.

SOURCE

## Source 1
source_surf                     = .false.        # source inside the medium, or source automatically moved exactly at the surface by the solver
xs                              = 96000.00      # source location x in meters
zs                              = 180000.00      # source location z in meters (zs is ignored if source_surf is set to true, it is replaced with the topography height)
## Source type parameters:
#  1 = elastic force or acoustic pressure
#  2 = moment tensor
# or Initial field type (when initialfield set in Par_file):
# For a plane wave including converted and reflected waves at the free surface:
#  1 = P wave,
#  2 = S wave,
#  3 = Rayleigh wave
# For a plane wave without converted nor reflected waves at the free surface, i.e. with the incident wave only:
#  4 = P wave,
#  5 = S wave
# For initial mode displacement:
#  6 = mode (2,3) of a rectangular membrane
source_type                     = 1
# Source time function:
# In the case of a source located in an acoustic medium,
# to get pressure for a Ricker in the seismograms, here we need to select a Gaussian for the potential Chi
# used as a source, rather than a Ricker, because pressure = - Chi_dot_dot.
# This is true both when USE_TRICK_FOR_BETTER_PRESSURE is set to .true. or to .false.
# Options:
#  1 = second derivative of a Gaussian (a.k.a. Ricker),
#  2 = first derivative of a Gaussian,
#  3 = Gaussian,
#  4 = Dirac,
#  5 = Heaviside (4 and 5 will produce noisy recordings because of frequencies above the mesh resolution limit),
#  6 = ocean acoustics type I,
#  7 = ocean acoustics type II,
#  8 = external source time function = 8 (source read from file),
#  9 = burst,
# 10 = Sinus source time function,
# 11 = Marmousi Ormsby wavelet
time_function_type              = 2
# If time_function_type == 8, enter below the custom source file to read (two columns file with time and amplitude) :
# (For the moment dt must be equal to the dt of the simulation. File name cannot exceed 150 characters)
# IMPORTANT: do NOT put quote signs around the file name, just put the file name itself otherwise the run will stop
name_of_source_file             = ""             # Only for option 8 : file containing the source wavelet
burst_band_width                = 0.             # Only for option 9 : band width of the burst
f0                              = 1.400e-02      # dominant source frequency (Hz) if not Dirac or Heaviside
tshift                          = 0.000e+00      # time shift when multi sources (if one source, must be zero)
## Force source
# angle of the source (for a force only); for a plane wave, this is the incidence angle; for moment tensor sources this is unused
anglesource                     = 0.00
## Moment tensor
# The components of a moment tensor source must be given in N.m, not in dyne.cm as in the DATA/CMTSOLUTION source file of the 3D version of the code.
Mxx                             = 1.000000       # Mxx component (for a moment tensor source only)
Mzz                             = -1.000000      # Mzz component (for a moment tensor source only)
Mxz                             = 0.000000       # Mxz component (for a moment tensor source only)
## Amplification (factor to amplify source time function)
factor                          = 1.000e+10      # amplification factor
## Moving source parameters
vx                              = 0.0            # Horizontal source velocity (m/s)
vz                              = 0.0            # Vertical source velocity (m/s)

interfaces.dat

# number of interfaces
 2
#
# for each interface below, we give the number of points and then x,z for each point
#
# interface number 1 (bottom of the mesh)
 2
 0.0 0.0
 480000.0 0.0
# interface number 5 (topography, top of the mesh)
 2
 0.0 180000.0
 480000.0 180000.0
#
# for each layer, we give the number of spectral elements in the vertical direction
#
# layer number 1
 40