daniel-koehn / denise-black-edition Goto Github PK

View Code? Open in Web Editor NEW

120.0 120.0 65.0 20.72 MB

2D time-domain isotropic (visco)elastic FD modeling and full waveform inversion (FWI) code for P/SV-waves

License: GNU General Public License v2.0

Makefile 0.12% C 23.18% Shell 0.07% Jupyter Notebook 74.66% Python 1.51% MATLAB 0.46%

denise-black-edition's People

Contributors

Stargazers

Watchers

denise-black-edition's Issues

Parsing bug in pyapi_denise.py

Hi everyone,

A bug in _parse_inp_file(self) function, starting from line 385 of pyapi_denise.py

The parameters file *.inp alongside with single argument per line contains multiple constants defined in the same line, e.g. GRAD1,GRAD2,GRAD3,GRAD4 = 15, 26, 19, 5. This case isn't taken into account in the parser and it attempts to create an argument of the Denise class with commas.

d = api.Denise()
# What we want to do
d.GRAD1 = 15
d.GRAD2 = 26
d.GRAD3 = 19
d.GRAD4 = 5
# What we actually have to do as workaround
d.__dict__["GRAD1,GRAD2,GRAD3,GRAD4"] = 15, 26, 19, 5

To fix this, need to tell first that there are many args per line, set each of them and attributes, and also be able to write them back to file into the correct location. Or make each parameter on its own line in the *.inp file.

Simulation is unstable in PSV mode

Hello, I have a question, regarding the work of FWI. In PSV mode I get error(Simulation is unstable !). Why do I get this, if CFL condition is satisfied and I have small enough dt?
Model on this iteration also looks reasonable. If you need more data from me, please ask.

There is a log file below.

Calculate residuals

iteste = 3 itests = 3 step1 = 1 step2 = 0 eps_scale = 3.375000e-02 countstep = 2 stepmax= 6 scalefac = 2.000000e+00 MYID = 0 L2t[1] = 1.148769e-16 L2t[2] = 1.112919e-16 L2t[3] = 1.081774e-16
MYID = 0 pimaxr = 4.228599e+03 gradmaxr = 3.681690e+05
MYID = 0 umaxr = 2.412584e+03 gradmaxr_u = 1.170411e+05
MYID = 0 rhomaxr = 2.555284e+03 gradmaxr_rho = 1.396093e+05

**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. 3 for shot 1 of 51 (rel. step length 0.03375000)
Message from function wavelet written by PE 0
1 source positions located in subdomain of PE 0
have been assigned with a source signal.

Time step: 2776; pvy: -nan
Message from PE 0
R U N - T I M E E R R O R:
Simulation is unstable !
...now exiting to system.
Message from PE 4
R U N - T I M E E R R O R:
Simulation is unstable !
...now exiting to system.

I tried to make my issue reproducible with data itself. Please download the archive and try to reproduce the issue.
https://drive.google.com/file/d/1fhIx2rlPzHXRfql-gPEr6HV9l9a3LY63/view?usp=sharing

Visco-Elastic modelling

Hi Dr. Daniel,

I am trying to model visco-elastic wave by changing Qs and Qp，but the modelling data do not change with Qs and Qp. I am sure "L=1" in my modelling. Why that happen? Looking forward to your reply. Thank You!!

Best
Liu

SegFault running example notebook and Python in general

Hello
I'm new to the world of geophysics, and was trying to use DENISE as my first FWI tool. When trying to run the example notebook, the forward modeling works just fine, but when doing FWI proper I get a segmentation fault error when calculating the gradients.

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

The error (part of it) is in the image attached. I get the same thing when running a simple model I created, no matter if from the command line or inside a notebook. Setting SWS_TAPER_GRAD_HOR=0 gives the same error in a later step, after

writing merged model file to /home/delta/bernardo/TestFWI_Denise/marmousi_tutorial/jacobian/gradient_Test_p_u.old Opening model files: /home/delta/bernardo/TestFWI_Denise/marmousi_tutorial/jacobian/gradient_Test_p_u.old.??? ... finished. Copying... ... finished. Use ximage n1=174 < /home/delta/bernardo/TestFWI_Denise/marmousi_tutorial/jacobian/gradient_Test_p_u.old label1=Y label2=X title=/home/delta/bernardo/TestFWI_Denise/marmousi_tutorial/jacobian/gradient_Test_p_u.old to visualize model.

Interestingly, running the steps listed in quickstart_MARMOUSI.txt works with no issues.

Any help is appreciated

Error: external source file can no be opened

Hi, Daniel,

The error in the first screen-shot comes up when I was trying to use an external source signal read from a file. The problem was solved by changing 'signal_wave' to 'SIGNAL_FILE' in Line 29, wavelet.c file as indicated in the second screen-shot. Can't figure out why.

Regards,

Ray

unable to run shot parallelization

When I try to use the shot parallelization feature using
d.fwi(model_init, src, rec, run_command = "mpirun -np 24"),
it doesn't run for the acoustic case.

Segmentation fault with READREC==2

Setting READREC==2 causes a segmentation fault.

When doing forward modeling, this occurs at

DENISE-Black-Edition/src/PSV/FD_PSV.c

Line 387 in 74b57f2

if (SEISMO) free_imatrix(acq.recpos,1,3,1,ntr_glob);

When doing FWI, this occurs at

DENISE-Black-Edition/src/PSV/FWI_PSV.c

Line 928 in 74b57f2

if (SEISMO) free_imatrix(acq.recpos,1,3,1,ntr_glob);

QUELLTYPB=4 fails to do FWI in AC problem

When I do SEISMO=2, QUELLTYPB=4, The LBFGS method fails to find the step length.

Bug in acoustic l-bfgs optimization

Currently, acoustic FWI (MODE=1, PHYSICS=2) does only converge, when using PCG optimization (GRAD_METHOD=1). The l-bfgs optimization implementation (GRAD_METHOD=2) does not converge.

I suspect there is a bug when using an even number of parameter classes, because the l-bfgs implementation works in GERMAINE for the mono-parameter Vp inversion and in DENISE with 3-parameter Vp, Vs, Density inversion.

Issues about boundary conditions.

@daniel-koehn
Hello, I am facing an issue while simulating an acoustic wave using a VSP geometry. I have observed that the CPML boundary conditions are not yielding the expected results. My model has a shape of 1280 (NX) * 610 (NY), with a DH of 5.0 and a PHYSICS parameter set to 2. To understand CPML, I have been reading research papers. During my experimentation, I tried several parameters for my geometry, but I discovered that only when FW=4, I get some positive results. However, there are still some issues that need to be resolved. I also attempted to use other parameters, such as npower and k_max_PML, but it appears that they have little effect on the results.
WaveField_Simulation.webm

Error opening field data in DENISE despite reading with Madagascar-based FWI and RTM

Hello, @daniel-koehn

First of all, I wanted to express my gratitude to you for recommending the NGP Course “Seismic Full Waveform Inversion” that has helped me get to this point.

I am using DENISE for FWI and RTM based on field seismic data. To prepare the data, I sorted my shot data and velocity field in Madagascar and SegyIO. Following this, I used Madagascar and Python to convert my segy data into su format. I normalized my shot data trace by trace using Madagascar and I ensured that the maximum amplitude was 1, as observed in the examples provided. After these preprocessing steps, I obtained the Source.dat and receiver.dat files and estimated the frequency of the source wavelet.

However, when I try to execute DENISE on my data, I receive an error message stating “Field data could not be opened”. I have checked the data and was able to read and plot it using Madagascar, so I guess that the issue is not with the data.

When I checked the su file with Madagascar, it appears as shown in the attached picture. It appears that both of these files have incorrect headers. The top file is the forward model by DENISE, and the bottom file is my field data.

Could you please provide some guidance on what might be causing this issue and how it could be resolved? Any help would be greatly appreciated. If you need any additional information regarding my data, please do not hesitate to let me know. I would be happy to provide more details. Alternatively, you can also email me at [email protected].

Best regards

Thank you!

Case INV_STF==1 is not seen in stf.c

Hi, Daniel,
I noticed that INV_STF==1 doesn't appear in your code of source wavelet inversion, stf.c. Does it mean that when INV_STF ==1, the time traces picked is controlled by the time damping?

Regards,

Ray

Save custom taper binaries to par/ directory

Hi everyone,

When creating custom taper geometries, the code searches for them in the directory from where the code was launched, all other data files are fetched from the respective par/ folder.

The experiment setup would be much more homogeneous if the data attributed to the same experiment would be stored in the same par/ folder.

Now you need to put tapers like this if starting simulations from this level (not from inside par folder):

bin/
...
par1/
par2/
taper.bin
taper_u.bin
taper_rho.bin

Would be better if like this (all experiment-specific tapers inside par folders):
Now you need to put tapers like this if starting simulations from this level (not from inside par folder):

bin/
...
par1/
----taper.bin
----taper_u.bin
----taper_rho.bin
par2/

Best regards,
Oleg

Creating rho,vp,vs files

Hello,

Hope you can help me.

I've created a velocity model to run on your tools and create shot gathers.

My problem is, I don't understand (and don't find in the manual) how to build a folder like: Marmousi-II/start folder.

Let's say that I have a 2D array of vp, do I just create a new file with ".vp" ending and put the values of the array inside it?

Another Q I have is, what is the difference between "start 1D" file and "margine" file?

Best,
Maayan

Seeking Guidance on Adapting DENISE Black Edition for Field Data Processing in Elastic Wave FWI and LSRTM

Hi Daniel,

I recently came across your excellent open-source code, DENISE Black Edition, which I believe will be extremely helpful for my master's thesis project, as I need to perform FWI and LSRTM of elastic waves.

Firstly, I would like to express my gratitude for your contribution to the open-source community by making your code available for use. I spent two weeks carefully reading the software manual that you wrote, and I tested the FWI and RTM examples you provided on OBC Marmousi. I also added my simple bit of code into your project (#42).

As part of my research, I need to work with field data for the FWI and LSRTM of elastic waves. I have been trying to read your program source code for the past two weeks, and I have carefully studied the application of actual data in your article "https://onlinelibrary.wiley.com/doi/10.1002/nsg.12097". However, I am having some difficulty understanding how to modify RTM-PSV.c, gra_obj_psv.c, stf_psv, and psv.c to incorporate real data into the framework.

With that said, I have a few questions for you:

Based on my preliminary understanding of your program, what do you think would be the main code changes required to use your code for actual data processing for the FWI and LSRTM of elastic waves?
Would it be possible for you to share some of your experience or code for elastic wave FWI and LSRTM for actual data processing?

Thank you very much for your time and expertise. I look forward to hearing back from you.

Best regards

Zhang

Full waveform inversion of 2-D acoustic wave

Hi Daniel,
First of all, thank you very much for sharing the program for everyone to use and learn; also thank you very much for your achievements in the professional field.
Recently, I was learning the code you shared. Your code is very practical. I encountered a problem when using the program to retrieve the full waveform of two-dimensional sound waves. The error report screenshot is as follows:

I used the P component of velocity as the observation data and set the density to a constant value: 1000kg / m3. When the inversion is carried out in four frequency ranges, the inversion results can be obtained in the first frequency range, and the above errors appear when the program continues to run. I looked at the current inversion results, the velocity value and density value are - Nan.
I want to ask you, have you ever encountered such a problem in the test, and how did you solve it?
Thank you again for your work and help. We look forward to your reply.

Multiple bugs in SH-FWI

The current implementation of SH-FWI in DENISE Black-Edition contains multiple bugs, which prevent a successfull FWI. Therefore, you have to use the stable SH-FWI code DENISE_SH

python IO

It would be great to have an easy to use interface to change input parameters from python 3.

Issues related to new gradient formulation GRAD_FORM=2

There are still some issues related to the gradients for the non-conservative stress-velocity formulation (Castellanos, 2014) - GRAD_FORM = 2:

The FWI only converges when using PCG optimization
The inversion of hydrophone data is currently not supported, because the source excitation for the backpropagtion of the data residuals requires some modifications
Compared to GRAD_FORM = 1 the ambiguity between Vs and density model is much higher. In case of the Marmousi-II model, structures from the Vs model dominate over density structures.

Combined FWI and gravity inversion

Steplength estimation during combined FWI and gravity inversion fails.

RUN_MULTIPLE_SHOTS issue

when RUN_MULTIPLE_SHOTS = 0 and Lines 223 in FD_PSV.c has been revised as "for (ishot=1;ishot<=1;ishot+=1)". The forward problem is correct, only shot 1 of 1 is simulated. But in inversion problem, shot 2~N of 1 are still be simulated.

QUELLTYPB=4 fails to do FWI in AC problem

Hi Dr. Daniel,

I tried to do FWI in AC problem. I first followed your guide 'QUICKSTART_Marmousi' and everything went well. Then I set output_of_seismograms_(SEISMO) = 2 and QUELLTYPB=4 so that only hydrophone data is used. This results in failing to estimate the steplength in less than one hour. Did I miss something?

I also found that the output time reversed residual seismograms are empty/0 by plotting.

Spatial 10th and 12th order FD operators

Currently, I only implemented the spatial 10th and 12th order FD operators for the elastic case.
The visco-elastic case requires a modification of update_s_visc_PML.c

Problem solved

Merge issues when no. of MPI-processes > 99

The extension of models written by each MPI-process

model.bin.%i%i

leads to merge problems when the number of MPI-processes is > 99.
Has to be replaced by

model.bin.%i.%i

Topography

Hey, dear Daniel,

I am a new user for your codes, I am wondering can I set a relief topography for the simulations?

Many thanks!

RTMOD not implemented

A module for the localization of microseismic events via reverse time modelling is not implemented in the current DENISE version.

problem with "make denise" code

Hello,
I am going through the QUICKSTART_Marmousi.txt and I have a problem with code "make denise"(step 3) . I tried all compiler options that exist in Makefile and it didn't work.
I working with Ubuntu os, gcc version 11.3.0, Open MPI 4.1.2, fftw version 3.3.8 .

currently I am using this configurations in Makefile:

On Desktop computer with LinuxMint 17, OpenMPI and gcc 4.8.2
CC=mpicc
LFLAGS=-lm -lcseife -lfftw3 -lstdc++
CFLAGS=-O3 -w -fno-stack-protector -D_FORTIFY_SOURCE=0
SFLAGS=-L./../libcseife
IFLAGS=-I./../libcseife -I./../include

this is file with the terminal output:
problem.txt

hope you can help.
thank you.

L5 (global correlation norm) misfit is negative. How to normalize it?

Hello, I am using the L5 or (global correlation) norm. As mentioned in the manual its implementation follows (Application of multi-source waveform inversion to marine streamer
data using the global correlation norm, Choi,Alkhalifah, 2012).

Originally it is negative.

I have a question: how to normalize it?

I attached the misfit file. And plots of the misfit functions on attached picture.

seis_fwi_log_L5.txt

New Gradient formulations GRAD_FORM

Different gradient formulations can be defined in denise.c

GRAD_FORM == 1 is based on the classical gradient derivation from the stress-displacement isotropic elastic equations of motion and reformulated for the stress-velocity code (e.g. Köhn et al. 2012). This approach has the disadvantage, that the data residuals have to be numerically integrated and are consequently very sensitive to noise.

EFWI with new gradient formulation GRAD_FORM == 2 for the stress-velocity formulation with compliance matrix (Vigh et al. 2014) is currently not converging, most likely due to an incorrect source implementation. The density update is also not implemented yet.

Gradient GRAD_FORM == 3 is based on the stress-velocity formulation with elastic tensor.
Actually, this problem is not self-adjoint and therefore the gradients should be wrong. However, the P-wave and S-wave velocity updates seem to converge correctly for the Marmousi-II test problem.
The density model shows a strong ambiguity and trade-offs with the S-wave velocity model.

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

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:

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.
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).
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).
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.
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

Inversion with L1 norm

Hi Daniel,

I tried doing an inversion with L1 norm and the Marmousi model. It didn't work.
Have you had experience with using the L1 norm for inversion?

Thanks and Regards,
Arnab

There is no such file

https://github.com/daniel-koehn/DENISE-Benchmark/raw/master/Marmousi-II/start/marmousi_II_marine.vp

Exception ("merge: can't read model file !") in mergemod.c

Sometimes, during my using of Denise PSV I get following error ("merge: can't read model file !") in mergemod.c.
What can be the reasons for this?
I am using 12 nodes 32 cpu each. NPROCX=4,NPROCY=4

**Message from mergemod (printed by PE 0):
PE 0 starts merge of 16 model files
writing merged model file to ./fwi/ws_fwi_3_strategy_51/Overthrust_true/fld/model/modelTest_vs_stage_1_it_10.bin
Opening model files: ./fwi/ws_fwi_3_strategy_51/Overthrust_true/fld/model/modelTest_vs_stage_1_it_10.bin.??? ... finished.
Copying... ... finished.
Use
ximage n1=384 < ./fwi/ws_fwi_3_strategy_51/Overthrust_true/fld/model/modelTest_vs_stage_1_it_10.bin label1=Y label2=X title=./fwi/ws_fwi_3_strategy_51/Overthrust_true/fld/model/modelTest_vs_stage_1_it_10.bin
to visualize model.

PE 0 is writing model to
./fwi/ws_fwi_3_strategy_51/Overthrust_true/fld/model/modelTest_rho_stage_1_it_10.bin.0.0

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

writing merged model file to ./fwi/ws_fwi_3_strategy_51/Overthrust_true/fld/model/modelTest_rho_stage_1_it_10.bin
Opening model files: ./fwi/ws_fwi_3_strategy_51/Overthrust_true/fld/model/modelTest_rho_stage_1_it_10.bin.??? Message from PE 0
R U N - T I M E E R R O R:
merge: can't read model file !
...now exiting to system.

-rw-r--r-- 1 plotnips k1404 0 May 19 22:17 modelTest_rho_stage_1_it_10.bin
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.0.0
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.0.1
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.0.2
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.0.3
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.1.0
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.1.1
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.1.2
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.1.3
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.2.0
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.2.1
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.2.2
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.2.3
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.3.0
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.3.1
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.3.2
-rw-r--r-- 1 plotnips k1404 90K May 19 22:17 modelTest_rho_stage_1_it_10.bin.3.3

2D time domain acoustic FWI

Hi Daniel,

Is there any implementation of 2D time domain acoustic FWI for totally freshman to dive into the field of FWI?

Thank you very much.

I can't find DENISE.inp in folder /DENISE/par

I've done the installation by reading manual_denise.pdf. And section 5.4 "Running the program" says I should have a DENISE.inp and FWI_workflow.inp in folder /par. However I only have DENISE_marm_OBC.inp, FWI_workflow_marmousi.inp and RTM_workflow_marmousi.inp three .inp files, the same in https://github.com/daniel-koehn/DENISE-Black-Edition/tree/master/par.

I've changed the filename DENISE_marm_OBC.inp to DENISE.inp and FWI_workflow_marmousi.inp to FWI_workflow.inp, and ran

mpirun -np 8 ../bin/denise DENISE.inp FWI_workflow.inp

It gave an RUN-TIME ERROR output below:

[iampora@localhost par]$ mpirun -np 8 ../bin/denise DENISE.inp FWI_workflow.inp
 *******************************************************************************
 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=0: Only forward modeling is applied.


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


In FD_PSV MYID = 0, COLOR =0, MYID_SHOT = 0 
 This is the log-file generated by PE 0 

Message from PE 0
R U N - T I M E  E R R O R: 
Number of processors specified in the parameter file 
 and at command line (NP) differ !
...now exiting to system.
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0
[iampora@localhost par]$ mpirun -np 8 ../bin/denise DENISE.inp FWI_workflow.inp > DENISE.out
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0

So where should I get those two .inp files and how should I run the example program?

QUELLTYPB = 4 not supported

The elastic FWI of hydrophone data (QUELLTYPB = 4) is currently not supported.

LBFGS

Hi Daniel,
I feel puzzled about the normlization and denormalization of gradients in LBFGS. Why the gradients multiply C_vp and C_rho in both normalization and denormalization? It will make the backtracking line search fail.

/* Normalization of the gradient   */
/* ------------------------------- */
for (i=1;i<=NX;i=i+IDX){
   for (j=1;j<=NY;j=j+IDY){
      waveconv[j][i] = C_vp * waveconv[j][i];
   }
}
  
for (i=1;i<=NX;i=i+IDX){
   for (j=1;j<=NY;j=j+IDY){
	  gradp[j][i] = waveconv[j][i];
   }
}

/* ===================================================================================================================================================== */
/* ===================================================== GRADIENT rho ================================================================================== */
/* ===================================================================================================================================================== */

/* Normalization of the gradient   */
/* ------------------------------- */
for (i=1;i<=NX;i=i+IDX){
   for (j=1;j<=NY;j=j+IDY){
      waveconv_rho[j][i] = C_rho * waveconv_rho[j][i];
   }
}

     /* Denormalize Gradients */
     for (i=1;i<=NX;i=i+IDX){
        for (j=1;j<=NY;j=j+IDY){
            
           waveconv[j][i] = waveconv[j][i] * C_vp;
	   waveconv_rho[j][i] = waveconv_rho[j][i] * C_rho;

        }
     }

Some comfuse about SH modelling

Hi Dr Kiehn,

When we modelling the SH wave in Denise-black-edition. How do we set the SH source in the source file？Should I select QUELLTYP=2? Look like we get the same modelling results using SH and PSV when we set QUELLTYP=2?

Thanks

DENISE-BE SH FWI currently only supports FD_ORDER = 2

Over the last 1.5 weeks, I tried to fix all of the bugs in the DENISE-BE SH module. DENISE-BE should now deliver the same results as the DENISE-SH code. I currently only checked the SH-FWI for the Ainos-SH field dataset: https://danielkoehnsite.wordpress.com/blog/ngp-course-seismic-full-waveform-inversion/ Due to time constraints, I could only implement the SH-FWI using FD_ORDER=2 operators. For higher order fd operators the SH-FWI will fail.

issue with using source wavelet in ASCII format

I have an issue with external wavelet file import (QUELLART=3) in DENISE MODE= 0 and = 1. Appears "R U N - T I M E E R R O R: Source file could no be opened !". File in ASCII format with one sample per line, sample interval and number of samples corresponds to the specified in the "DENISE***.inp" file. Log files attached.

DENISE.out.txt
US.log.1.txt