Zhaolun Liu
King Abdullah University of Science and Technology
This lab is for implementation of the 3D topographic wave equation dispersion inversion (3D TWD). It includes
Even though the above package number from 1 to 3 are open source packages, it may be not suitable for 3DTWD or may include some bugs. So, some changes are needed for these packages. For example, I modified the "write_output_SU" in SPECFEM3D to scale the coordinates of receivers by 1000 to get a more accurate value. I modified a lot in SeisFlows, and the main changes are listed as follows,
A: You can open "3dfoothill/output.optim" to see if the misfit is decreasing. Also, you can check the synthetic data, dispersion image (by Radon transform+fft) and picked dispersion curve in "3dfoothill/output/traces_000*".
obs/000000/dispObs.r: picked dispersion curve
obs/000000/fk_dataobs.r: dispersion image
obs/000000/obs.r: observed data
syn/000000/dispPre.r: picked dispersion curve
syn/000000/fk_data.r: dispersion image
syn/000000/syn.r: predicted data
The above files have both a header and a binary file.
The inverted model is in "3dfoothill/output/model_000*". You can check the model by generating the VTK file from the inverted model using "xcombine_vol_data_vtk" which should in the "SPECFEM3D/bin" For example: you can use the attached script.
./create_vs_vtk.sh vs 21 input_folder
which requires:
1, you have "bin/xcombine_vol_data_vtk", "DATA/Par_file" in you folder.
2, copy "OUTPUT_FILES/DATABASES_MPI/proc00000*_external_mesh.bin" to your local folder.
A: The model is decomposed into 22 parts so that there are 22 files for Vs. The dimension information is stored in "proc0000**_external_mesh.bin", which can be found in "scratch/solver/000000/OUTPUT_FILES/DATABASES_MPI/". For example you can use an attached matlab script modelchange.m ``modelchange.m'' to read and modify them.
You can check the model by generating the VTK file from the inverted model using xcombine_vol_data_vtk which should in the "SPECFEM3D/bin" For example: you can use the attached script.
./create_vs_vtk.sh vs 21 input_folder
which requires: 1, you have "bin/xcombine_vol_data_vtk", "DATA/Par_file" in you folder. 2, copy "OUTPUT_FILES/DATABASES_MPI/proc00000*_external_mesh.bin" to your local folder.
A: I encountered this phenomenon before. I suppose that the inversion has reached a local minimum. But you can increase the trial number by setting STEPMAX in line 58 in parameters.py (shown in Fig. 2).
A: The test is used just one shot so that the result is not so good. You can use all the 80 sources by setting NSRC=80 in line 19 of parameters.py. If you increase NSRC, please also set NTASK in line 74, parameters.py.
If you want to do forward modeling by yourself. You can set DATA=' ' in paths.py. The code will automatically check if DATA is empty. If DATA is empty, it will do forward modeling from the true model. Otherwise, you can set it to "path/obsdata_3000"
450.000000 380.000000
450.000000 1330.000000
450.000000 2280.000000
450.000000 3230.000000
1400.000000 380.000000
1400.000000 1330.000000
1400.000000 2280.000000
1400.000000 3230.000000
2350.000000 380.000000
2350.000000 1330.000000
2350.000000 2280.000000
2350.000000 3230.000000
3300.000000 380.000000
3300.000000 1330.000000
3300.000000 2280.000000
3300.000000 3230.000000
4250.000000 380.000000
4250.000000 1330.000000
4250.000000 2280.000000
4250.000000 3230.000000
5200.000000 380.000000
5200.000000 1330.000000
5200.000000 2280.000000
5200.000000 3230.000000
6150.000000 380.000000
6150.000000 1330.000000
6150.000000 2280.000000
6150.000000 3230.000000
I suppose that these are the locations of the sources, but there are only 28 of them, not 80.
A: The "source.dat" file is not used in the inversion. FORCESOLUTION_0000** is actually used in the inversion.
A: For a new dataset: you need to do the followings:
a) define the mesh by using CUBIT. An example can be found in "specfem3d/EXAMPLES/Mount_StHelens/"
b) generate a homogeneous gll model file:
You can do a forward modeling test by using the model parameters described by mesh properties. For example, set "MODEL=default" in the Par_file. Then the gll model will be generated in "OUTPUT_FILES/DATABASES_MPI/proc00000**_vp/vs/rho/external_mesh.bin"
c) define the initial model
You can define a new model based the above gll model, which can be done by a matlab script that I sent you in the previous email.
d) define STATIONS and FORCESOLUTION
e) define the parameters for calculating the dispersion curve as shown
in Fig. 3.
f) VSMAX VSMIN in the parameters.py is the clip for s-veloicty model
For field data, "lzdispertest" need to be modified so that it can read the manually-picked dispersion curves, not calculate the dispersion curve within the code.
A: If you want to do modeling, you may try to change "WORKFLOW='test_forward' " in parameters.py.
A: The documents of SPECFEM3D shows how to define the location of the stations. Each line represents one station in the following format:
Station Network Latitude(degrees) Longitude(degrees) Elevation(m) burial(m)
The solver xspecfem3D filters the list of stations in file DATA/STATIONS to exclude stations that are not located within the region given in the Par_file (between LATITUDE_MIN and LATITUDE_MAX and between LONGITUDE_MIN and LONGITUDE_MAX). The filtered file is called DATA/STATIONS_FILTERED. Elevation and burial are generally applicable to geographical regions. Burial is measured down from the top surface. For other problems in other fields (ultrasonic testing, medical imaging etc...), it may be confusing. We generally follow either of the following procedures for those kind of problems: Procedure 1:
- Always put the origin on the top of the model.
- Let’s say you want to place two receivers at (x1,y1,z1) and (x2,y2,z2). Your STATIONS file should look like:
BONE GR y1 x1 0.00 -z1
BONE GR y2 x2 0.00 -z2
Depending on the geometry of the models it may not work in some cases.
For my example, the burial is zeros so that you find the depth is zero. For more information, you can refer to Chapter 5 of User Manual of SPECFEM 3D Cartesian
A: It means you only used 20 shots by using SYSTEM="slurm_lg".
NSRC=80 # number of sources
NTASK=80 # number of tasks or sources
NODESIZE=36
SLURMARGS='-p BDW'
NPROC=22 # number of processers model decomposition
In the set-up above, I do not know how many cores are used.
If I kill one of the many jobs in the middle of inversion process, will the SeisFlow inversion continue by using the remaining cores?
A:
The current job will use 80*22 cores. The following command how the parameter are used:
sbatch -A k1046 -job-name=3dfoothill
-nodes=math.ceil(PAR.NPROC/float(PAR.NODESIZE)) -ntasks-per-node=PAR.NODESIZE -ntasks=PAR.NPROC -time=PAR.TASKTIME -array=0-PAR.NTASK-1 -output /lustre/project/k1046/liuz0b/newspecfem/test/3dfoothill/output.slurm/%A_%a
/project/k1046/liuz0b/newspecfem/seisflows-master/seisflows/system//wrappers/run
/lustre/project/k1046/liuz0b/newspecfem/test/3dfoothill/output solver setup
which comes from "system/slurm.lg.py"
It allocate PAR.NTASK*
(math.ceil(PAR.NPROC/float(PAR.NODESIZE))=1 in this case) nodes and each node use PAR.NPROC core.
By using SYSTEM="slurm_lg", it cannot implement the job submission that you described. for the job submission like what you said, I don't know how to do it in seisflow. Maybe, you can see other scripts in "seisflow/system/*.py". It has "slurm_sm" but I don't think that it can implement what you want.
Because it use the "array" in sbatch. You should kill both the main job and the array. The job ID name for the array usually looks like jobID_0, jobID_1, ..., jobID_NTASK-1. you can kill it by the command "scancel jobID". If you just kill jobID_0, the remaining job in the array will continue to run.
The NTASK actually control the number of source, not NSRC.