Seeking Advice on FWI using DENISE with OBN Field Data: Challenges and Simulation Instability

[zhangxiaoshuotttt]

Hello everyone @daniel-koehn @ovcharenkoo @ADharaUTEXAS123007 @pplotn, I hope you’re doing well. I would like to seek your advice and assistance regarding Full Waveform Inversion (FWI) using DENISE with OBN field data. I have learned that you all are familiar with and have expertise in DENISE and FWI, which is why I am reaching out for your guidance. I have encountered some challenges and would greatly appreciate any assistance you can provide. Thank you in advance for your help!

Here is a brief overview of what I have accomplished in the past month:

1. First, I selected a source line that closely aligns with the receiver line, ensuring that the receiver interval remains constant. Then, I projected the receiver line onto the source line.

2. I performed a 3-component pre-stack time migration using the OBN data. During this process, I obtained Vp and Vs values in the time domain, which I subsequently converted into the depth domain. Here is the overview of my acquisition (to save time and for testing purposes, I used only 16 shots).

3. I utilized the raw data, meaning that I did not apply any conventional seismic preprocessing to my data (though I’m uncertain if this is the correct approach).

4. To generate synthetic data, I used the velocity field and density (rho) field. I’m attaching both my synthetic data and the actual field data for your reference. I observed that the Marmousi example also includes a water layer. However, when I set the water layer in the Vs field, my water depth is approximately 80 meters (same as the receiver positions), while my source position is at around 6 meters. Unfortunately, I couldn’t obtain the desired synthetic data until I moved the source position to 12.5 meters.
   
   
   

5. During my attempts to achieve FWI using DENISE, the first epoch performed well. However, upon entering the second epoch, I encountered an error stating “Simulation is unstable.” I understand that the reason behind this issue is likely due to my lack of experience with FWI, specifically in obtaining the gradient field. As a newcomer to this field, I am struggling to analyze the exact cause of this occurrence.

 *******************************************************************************
 This is program DENISE Black-Edition                                           
 Parallel 2-D elastic Finite Difference FWI code                                
                                                                                
 Forward/FWI/RTM codes written by D. Koehn and D. De Nil                        
 2D isotropic PSV forward code partly based on FDVEPS written by T. Bohlen      
 Institute of Geosciences, Kiel University, Germany                           

 See README.md file and LICENSE.md for redistribution conditions.               
 *******************************************************************************


 **Message from check_model_phys (printed by PE 0):

 -----------------------  DENISE operation mode  ----------------------
 MODE=1: Time-domain FWI is applied.


 -----------------------  DENISE Physics  ----------------------
 PHYSICS=1: Solve 2D isotropic elastic PSV problem.


 This is the log-file generated by PE 0 


 **Message from initprocs (printed by PE 0):
 Size of subarrays in gridpoints:
 IENDX= 628
 IENDY (vertical) = 250

 **Message from initprocs (written by PEEEP 0):
 Processor locations in the 2D logical processor array
 MYID 	 POS(1):left,right 	 POS(2): top, bottom
 0 		 0: 2,1 		 0: 9,3 

 **Message from write_par (printed by PE 0):

------------------------- Processors ------------------------
 Number of PEs in x-direction (NPROCX): 3
 Number of PEs in vertical direction (NPROCY): 4
 Total number of PEs in use: 12

 ----------------------- Discretization  ---------------------
 Number of gridpoints in x-direction (NX): 1884
 Number of gridpoints in y-direction (NY): 1000
 Grid-spacing (DH): 6.250000e+00 meter
 Time of wave propagation (T): 4.000000e+00 seconds
 Timestep (DT): 1.000000e-03 seconds
 Number of timesteps: 4000 

 ------------------------- FD ORDER -----------------------------
 FDORDER = 8
 MAXRELERROR = 1

 ------------------------- SOURCE -----------------------------
 reading source positions, time delay, centre frequency 
 and initial amplitude from ASCII-file 
	./source/OBNtest.dat

 wavelet of source: 1st derivative of Gaussian 


 ------------------------- RECEIVER  --------------------------
 reading receiver positions from single file 
	./receiver/OBN_ReceiverFile.dat

 reference_point_for_receiver_coordinate_system:
 x=0.000000 	y=0.000000	 z=0.000000
 ------------------------- FREE SURFACE ------------------------
 free surface at the top of the model ! 

 ------------------------- CPML ---------------------
 width of absorbing frame is 10 gridpoints.
 CPML damping applied. 
 Damping velocity in the PML frame in m/s: 1500.000000 .
 Frequency within the PML frame in Hz: 20.000000 
 npower: 4.000000 
 k_max: 1.000000 
 No periodic boundary condition.
 ------------------------- MODEL-FILES -------------------------
 names of model-files: 
	 shear wave velocities:
	 start/OBN_model.vs
	 tau for shear waves:
	 start/OBN_model.ts
	 density:
	 start/OBN_model.rho
	 compressional wave velocities:
	 start/OBN_model.vp
	 tau for P-waves:
	 start/OBN_model.tp

 ------------------------- Q-APROXIMATION --------------------
 Number of relaxation mechanisms (L): 0
 The L relaxation frequencies are at:  
 Hz
 Value for tau is : 0.000010

 -----------------------  SEISMOGRAMS  ----------------------
 seismograms of  x- and y-component of particle velocity.
 output-files: 
 	su/OBN_RealData_x.su
	su/OBN_RealData_y.su
 The data is written in IEEE SU-format . 
 samplingrate of seismic data: 0.001000 s
 Number of samples per trace: 4000 
 ----------------------------------------------------------



 -----------------------  DENISE elastic specific parameters  ----------------------

 Maximum number of iterations: 100
 location of the measured seismograms : 
 	su/OBNTest/OBN_FinalRealdata

 INVMAT1=1: Inversion parameters are vp, vs and rho.
 QUELLTYPB=1: Inversion of x and y component.

 Shots used for step length estimation:
	 TESTSHOT_START = 1 
	 TESTSHOT_END = 31 
	 TESTSHOT_INCR = 10 

 Cosine Taper used : 
 	0
 Log file for misfit in each iteration step: 
 	LOG_TEST.dat 

 Output of inverted models to: 
 	model/DFB/modelTest 

 --------------- Gradient tapering -------------------
 SWS_TAPER_GRAD_VERT=1: Vertical taper applied.
 (GRADT1=21, GRADT2=25, GRADT3=490, GRADT4=500)

 SWS_TAPER_GRAD_HOR=1: Horizontal taper applied.
 (GRADT1=21, GRADT2=25, GRADT3=490, GRADT4=500, EXP_TAPER_GRAD_HOR=2.000000)

 SWS_TAPER_GRAD_SOURCES=0: No taper around the sources applied.

 SWS_TAPER_CIRCULAR_PER_SHOT=0: No taper around the sources applied.



 --------------- Gradient smoothing with 2D-Gaussian filter -------------------
 GRAD_FILTER=0: Gradients are not filtered.


 --------------- Limits of model parameters -------------------
 VPLOWERLIM = 1400.000000 		 VPUPPERLIM = 5000.000000 
 VSLOWERLIM = 0.000000 		 VSUPPERLIM = 4000.000000 
 RHOLOWERLIM = 0.000000 	 RHOUPPERLIM = 3000.000000 


 --------------- Optimization method -------------------
 GRAD_METHOD=2: LBFGS
 NLBFGS=20 


 --------------- Model smoothing -------------------
 MODEL_FILTER=0: vp and vs models are not filtered after each iteration step.


 --------------- Trace kill -------------------
 TRKILL=0: No trace kill is applied 


 --------------- Trace normalization -------------------
 NORMALIZE=0: No normalization of measured and synthetic seismograms.


 --------------- Reduce size of inversion grid -------------------
 Every 3 time sample is used for the calculation of the gradients.



 --------------- Step length estimation -------------------
 EPS_SCALE = 0.020000
 STEPMAX = 6
 SCALEFAC = 2.000000


 --------------- Reverse Time Modelling -------------------
 RTMOD=0: No Reverse Time Modelling applied.


 --------------- Gravity Modelling and Inversion -------------------
 No Gravity Modelling and Inversion applied. 	 GRAVITY=0 
 Boundary in x-direction [gridpoints] 			 NGRAVB = 500 
 Boundary in z-direction [m] 				 NZGRAV = 200000 
 Modelling and Inversion of Gravity Data. 		 GRAV_TYPE=1 
 Self-defined Model is used as Background Density. 	 BACK_DENSITY=2 
	 background density file:
	 gravity/background_density.rho 

 **************************************************************


 Reading receiver positions from single file: 
	./receiver/OBN_ReceiverFile.dat
 Message from function receiver (written by PE 0):
 Number of receiver positions found: 225

 **Message from main (printed by PE 0):
 Size of local grids: NX=628 	 NY=250
 Each process is now trying to allocate memory for:
 Dynamic variables: 		   3.13 MB
 Static variables: 		   3.76 MB
 Seismograms: 			   2.26 MB
 Buffer arrays for grid exchange:   0.07 MB
 Storage of forward modelled wavefields: 3994.71 MB
 Gradients, updated material parameters:  13.77 MB
 Residual seismograms: 	   3.39 MB
 Full seismograms: 	  10.30 MB
 Network Buffer for MPI_Bsend: 	   0.13 MB
 ------------------------------------------------ 
 Total memory required: 	 4021.22 MB.

 ... memory allocation for PE 0 was successfull.


 Reading source positions, time-shift, centre frequency 
 and amplitude from file: ./source/OBNtest.dat
 Number of source positions specified in ./source/OBNtest.dat : 16 
 Maximum frequency defined in ./source/OBNtest.dat: 2.00e+01 Hz
 All sources will be modelled individually because of RUN_MULTIPLE_SHOTS=1!
 Message from function sources (written by PE 0):

...reading model information from modell-files...
	 Vp:
	 start/OBN_model.vp

	 Vs:
	 start/OBN_model.vs

	 Density:
	 start/OBN_model.rho



 PE 0 is writing model to 
	start/OBN_model.denise.pi.0.0


 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  start/OBN_model.denise.pi 
 Opening model files: start/OBN_model.denise.pi.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < start/OBN_model.denise.pi  label1=Y label2=X title=start/OBN_model.denise.pi 
 to visualize model. 


 PE 0 is writing model to 
	start/OBN_model.denise.mu.0.0


 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  start/OBN_model.denise.mu 
 Opening model files: start/OBN_model.denise.mu.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < start/OBN_model.denise.mu  label1=Y label2=X title=start/OBN_model.denise.mu 
 to visualize model. 


 PE 0 is writing model to 
	start/OBN_model.denise.rho.0.0


 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  start/OBN_model.denise.rho 
 Opening model files: start/OBN_model.denise.rho.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < start/OBN_model.denise.rho  label1=Y label2=X title=start/OBN_model.denise.rho 
 to visualize model. 



 **Message from checkfd (printed by PE 0):
 Minimum and maximum P-wave and S-wave velocities within subvolumes: 
  MYID	 Vp_min(f=fc) 	 Vp_max(f=inf) 	 Vs_min(f=fc) 	 Vsmax(f=inf) 
 0 	 1.459000e+03 	 2.090969e+03 	 0.000000e+00 	 1.344821e+03 
 Global values for entire model: 
 Vp_max= 3.251180e+03 m/s 	 Vs_min=0.000000e+00 m/s 



 ------------------ CHECK FOR GRID DISPERSION --------------------
 To satisfactorily limit grid dispersion the number of gridpoints 
 per minimum wavelength (of S-waves) should be 6 (better more).
 Here the minimum wavelength is assumed to be minimum model phase velocity 
 (of S-waves) at maximum frequency of the source
 devided by maximum frequency of the source.
 Maximum frequency of the source is approximately    40.00 Hz
 The minimum wavelength (of S-waves) in the following simulation will
 be 0.000000e+00 meter.
 Thus, the recommended value for DH is 0.000000e+00 meter.
 You have specified DH= 6.250000e+00 meter.

W A R N I N G   M E S S A G E: 
 Grid dispersion will influence wave propagation, choose smaller grid spacing (DH).
 

 ----------------------- CHECK FOR STABILITY ---------------------
 The following simulation is stable provided that

 	 p=cmax*DT/DH < 1/(sqrt(2)*gamma),

 where cmax is the maximum phase velocity at infinite frequency
 and gamma = sum(|FD coeff.|)
 In the current simulation cmax is  3251.18 m/s .

 DT is the timestep and DH is the grid size.

 In this simulation the stability limit for timestep DT is 1.009956e-03 seconds .
 You have specified DT= 1.000000e-03 s.
 The simulation will be stable.


 ----------------------- ABSORBING BOUNDARY ------------------------
 Width (FW) of absorbing frame should be at least 10 gridpoints.
 You have specified a width of 10 gridpoints.
=========================================== 
       FWI-stage 1 of 4 
=========================================== 
 Density is inverted from iteration step 0.

 Vp is inverted from iteration step 0.

 Vs is inverted from iteration step 0.

 Qs is inverted from iteration step 0.


 Smoothing (spatial filtering) of the gradients: 
  	SPATFILTER=0: Gradients are not smoothed.


 --------------- Frequency filtering -------------------
 TIME_FILT=1: Time domain filtering is applied 
 Applying FWI with lowpass filter for corner frequency 2.000000 Hz
 Order of lowpass filter is:	6



 --------------- Time windowing and damping -------------------
 TIMEWIN=0: No time windowing and damping is applied 


 --------------- termination of the program -------------------
 Misfit change during the last two iterations is smaller than 1.000000 percent.



 --------------- Energy preconditioning -------------------
EPRECOND = 3 - Hessian approximation by Plessix & Mulder (2004) 

EPRECOND = 3 - energy preconditioning deactivated 



 --------------- Gradient calculation -------------------
 Misfit function:
   LNORM==1 corresponds to L1 Norm
   LNORM==2 corresponds to L2 Norm
   LNORM==3 corresponds to Cauchy
   LNORM==4 corresponds to SECH
   LNORM==5 corresponds to global correlation
   LNORM==6 corresponds to envelope objective function

 Used LNORM=2
 N_ORDER=0





 --------------- No STF inversion  -------------------



 --------------- TRACE normalization  -------------------

 NORMALIZE=1



 --------------- No Offset mute  -------------------

 --------------- Scale density update  -------------------
 SCALERHO = 0.500000

 --------------- Scale Qs update  -------------------
 SCALEQS = 1.000000

---------------- No Joint Inversion ----------------


                   TDFWI ITERATION 1 	 of 100 


 **Message from matcopy (written by PE 0):
 Copy material properties at inner boundaries ... 
 finished (real time: 0.00 s).
Vp_avg = 2486 	 Vs_avg = 1473 	 rho_avg = 2098 
  Message from function wavelet written by PE 0 
 1 source positions located in subdomain of PE 0 
 have been assigned with a source signal. 

 PE 0 outputs source time function in SU format to start/OBN_model_source_signal.0.su.shot1 
 
 *****  Starting simulation (forward model) for shot 1 of 16  ********** 
-------------------  
Calculate residuals  
-------------------  
 PE 0 is writing 225 seismograms (vx) to
	 su/OBN_RealData_x.su.shot1.it1 
 PE 0 is writing 225 seismograms (vy) to
	 su/OBN_RealData_y.su.shot1.it1 

 *****  Starting simulation (adjoint wavefield)  ********** 
 Message from function wavelet written by PE 0 
 1 source positions located in subdomain of PE 0 
 have been assigned with a source signal. 

 PE 0 outputs source time function in SU format to start/OBN_model_source_signal.0.su.shot2 
 
 *****  Starting simulation (forward model) for shot 2 of 16  ********** 
-------------------  
Calculate residuals  
-------------------  
.
.
.
.
 PE 2 outputs source time function in SU format to start/OBN_model_source_signal.2.su.shot16 
 
 *****  Starting simulation (forward model) for shot 16 of 16  ********** 
-------------------  
Calculate residuals  
-------------------  

 *****  Starting simulation (adjoint wavefield)  ********** 

 **Message from taper_grid (printed by PE 0):
 Coefficients for gradient taper are now calculated.


 PE 0 is writing model to 
	taper/taper_coeff_vert.bin.0.0


 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  taper/taper_coeff_vert.bin 
 Opening model files: taper/taper_coeff_vert.bin.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < taper/taper_coeff_vert.bin  label1=Y label2=X title=taper/taper_coeff_vert.bin 
 to visualize model. 

 **Message from taper_grid (printed by PE 0):
 Coefficients for gradient taper are now calculated.
21 	 25 	 490 	 500 
4 	 10 	 5 


 PE 0 is writing model to 
	taper/taper_coeff_hor.bin.0.0


 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  taper/taper_coeff_hor.bin 
 Opening model files: taper/taper_coeff_hor.bin.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < taper/taper_coeff_hor.bin  label1=Y label2=X title=taper/taper_coeff_hor.bin 
 to visualize model. 

 **Message from taper_grid (printed by PE 0):
 Coefficients for gradient taper are now calculated.


 PE 0 is writing model to 
	taper/taper_coeff_vert.bin.0.0


 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  taper/taper_coeff_vert.bin 
 Opening model files: taper/taper_coeff_vert.bin.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < taper/taper_coeff_vert.bin  label1=Y label2=X title=taper/taper_coeff_vert.bin 
 to visualize model. 

 **Message from taper_grid (printed by PE 0):
 Coefficients for gradient taper are now calculated.
21 	 25 	 490 	 500 
4 	 10 	 5 


 PE 0 is writing model to 
	taper/taper_coeff_hor.bin.0.0


 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  taper/taper_coeff_hor.bin 
 Opening model files: taper/taper_coeff_hor.bin.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < taper/taper_coeff_hor.bin  label1=Y label2=X title=taper/taper_coeff_hor.bin 
 to visualize model. 

 **Message from taper_grid (printed by PE 0):
 Coefficients for gradient taper are now calculated.


 PE 0 is writing model to 
	taper/taper_coeff_vert.bin.0.0


 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  taper/taper_coeff_vert.bin 
 Opening model files: taper/taper_coeff_vert.bin.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < taper/taper_coeff_vert.bin  label1=Y label2=X title=taper/taper_coeff_vert.bin 
 to visualize model. 

 **Message from taper_grid (printed by PE 0):
 Coefficients for gradient taper are now calculated.
21 	 25 	 490 	 500 
4 	 10 	 5 


 PE 0 is writing model to 
	taper/taper_coeff_hor.bin.0.0


 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  taper/taper_coeff_hor.bin 
 Opening model files: taper/taper_coeff_hor.bin.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < taper/taper_coeff_hor.bin  label1=Y label2=X title=taper/taper_coeff_hor.bin 
 to visualize model. 

 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  jacobian/OBN/jacobian_Test_p_vp.old 
 Opening model files: jacobian/OBN/jacobian_Test_p_vp.old.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < jacobian/OBN/jacobian_Test_p_vp.old  label1=Y label2=X title=jacobian/OBN/jacobian_Test_p_vp.old 
 to visualize model. 

 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  jacobian/OBN/jacobian_Test_p.old 
 Opening model files: jacobian/OBN/jacobian_Test_p.old.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < jacobian/OBN/jacobian_Test_p.old  label1=Y label2=X title=jacobian/OBN/jacobian_Test_p.old 
 to visualize model. 

 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  jacobian/OBN/jacobian_Test_c.old 
 Opening model files: jacobian/OBN/jacobian_Test_c.old.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < jacobian/OBN/jacobian_Test_c.old  label1=Y label2=X title=jacobian/OBN/jacobian_Test_c.old 
 to visualize model. 

 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  jacobian/OBN/jacobian_Test_p_vs.old 
 Opening model files: jacobian/OBN/jacobian_Test_p_vs.old.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < jacobian/OBN/jacobian_Test_p_vs.old  label1=Y label2=X title=jacobian/OBN/jacobian_Test_p_vs.old 
 to visualize model. 

 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  jacobian/OBN/jacobian_Test_p_u.old 
 Opening model files: jacobian/OBN/jacobian_Test_p_u.old.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < jacobian/OBN/jacobian_Test_p_u.old  label1=Y label2=X title=jacobian/OBN/jacobian_Test_p_u.old 
 to visualize model. 

 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  jacobian/OBN/jacobian_Test_c_u.old 
 Opening model files: jacobian/OBN/jacobian_Test_c_u.old.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < jacobian/OBN/jacobian_Test_c_u.old  label1=Y label2=X title=jacobian/OBN/jacobian_Test_c_u.old 
 to visualize model. 

 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  jacobian/OBN/jacobian_Test_p_mrho.old 
 Opening model files: jacobian/OBN/jacobian_Test_p_mrho.old.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < jacobian/OBN/jacobian_Test_p_mrho.old  label1=Y label2=X title=jacobian/OBN/jacobian_Test_p_mrho.old 
 to visualize model. 

 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  jacobian/OBN/jacobian_Test_p_rho.old 
 Opening model files: jacobian/OBN/jacobian_Test_p_rho.old.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < jacobian/OBN/jacobian_Test_p_rho.old  label1=Y label2=X title=jacobian/OBN/jacobian_Test_p_rho.old 
 to visualize model. 

 **Message from mergemod (printed by PE 0):
 PE 0 starts merge of 12 model files 

 writing merged model file to  jacobian/OBN/jacobian_Test_c_rho.old 
 Opening model files: jacobian/OBN/jacobian_Test_c_rho.old.???  ... finished. 
 Copying... ... finished. 
 Use 
 ximage n1=1000 < jacobian/OBN/jacobian_Test_c_rho.old  label1=Y label2=X title=jacobian/OBN/jacobian_Test_c_rho.old 
 to visualize model. 
MYID = 0 	 pimaxr = 3.251180e+03 	 gradmaxr = 0.000000e+00 
MYID = 0 	 umaxr = 2.038207e+03 	 gradmaxr_u = 0.000000e+00 
MYID = 0 	 rhomaxr = 2.265330e+03 	 gradmaxr_rho = 0.000000e+00 


 **Message from matcopy (written by PE 0):
 Copy material properties at inner boundaries ... 
 finished (real time: 0.00 s).

 *****  Starting simulation (test-forward model) no. 2 for shot 1 of 16 (rel. step length 0.02000000) 
 Message from function wavelet written by PE 0 
 1 source positions located in subdomain of PE 0 
 have been assigned with a source signal. 
Message from PE 6
R U N - T I M E  E R R O R: 
 Simulation is unstable !
...now exiting to system.
Message from PE 7
R U N - T I M E  E R R O R: 
 Simulation is unstable !
...now exiting to system.
Message from PE 8
Abort(1) on node 6 (rank 6 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 6
Abort(1) on node 7 (rank 7 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 7

I have updated my files to include detailed settings. Could you please provide me with some advice based on this information?Once again, I would like to express my deepest gratitude for your valuable time, expertise, and assistance. Your guidance will be immensely valuable to overcome the challenges I’m facing. Thank you all for your kind support.
DENISE_OBN_fIeldData_Flow.txt
FWI_workflow_marmousi.txt

Best regards,

Zhang

[zhangxiaoshuotttt]

I believe I have identified the main cause of this problem, which is related to the taper parameter. By adjusting the taper parameter, I was able to successfully resolve this issue.