daniel-koehn / denise-black-edition Goto Github PK
View Code? Open in Web Editor NEW2D time-domain isotropic (visco)elastic FD modeling and full waveform inversion (FWI) code for P/SV-waves
License: GNU General Public License v2.0
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
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
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.
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
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
When doing FWI, this occurs at
DENISE-Black-Edition/src/PSV/FWI_PSV.c
Line 928 in 74b57f2
Steplength estimation during combined FWI and gravity inversion fails.
FP
is declared at:
DENISE-Black-Edition/src/PSV/obj_psv.c
Line 24 in 74b57f2
and then used uninitialized when READREC==2
at:
DENISE-Black-Edition/src/PSV/obj_psv.c
Line 57 in 74b57f2
A module for the localization of microseismic events via reverse time modelling is not implemented in the current DENISE version.
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!
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
Li
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
Hello Dr. Daniel,
My name is Raul Leon. I have been using DENISE for some synthetic acoustic inversion problems and now I am trying to run an acoustic FWI using real data. Since DENISE creates the seismograms using an old version of seismic unix, I tried to run the inversion using ascii data-format but I got the following error:
Message from PE 1
R U N - T I M E E R R O R:
Field data file could not be opened !
...now exiting to system.
Message from PE 0
R U N - T I M E E R R O R:
Field data file could not be opened !
...now exiting to system.
Based on that, it seems that FWI workflow using ascii data-format is not available yet. So, in order to run the inversion using the seismic unix data format, which version of SU would you recommend me to use?
PE 0 is writing 648 seismograms of pressure to
su/DENISE_L0970_Tiber_p.su.shot1.it1
==================================================================================
***** Starting simulation (adjoint wavefield) **********
==================================================================================
Time step: 121; 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.
From the output file I know that my simulation should be stable so that makes me believe the error could be related to the seismic unix version that I am using to convert from segy to su. Could you give me any recommendation to avoid these types of problems?
I really appreciate your help.
Hello daniel, is there any existing parameter in the code that can set the lower surface as a free boundary? If not, how should I modify the code file?
It would be great to have an easy to use interface to change input parameters from python 3.
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 are still some issues related to the gradients for the non-conservative stress-velocity formulation (Castellanos, 2014) - GRAD_FORM = 2:
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
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.
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.
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.
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
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.
When I do SEISMO=2, QUELLTYPB=4, The LBFGS method fails to find the step length.
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.
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
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
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!
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
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
Dear Daniel,
I hope this email finds you well. My name is Chen. I am writing to inquire about the current version of the program and whether it supports SH RTM . I attempted to run the program recently, but it seems that I encountered difficulties in getting it to function properly for SH RTM purposes.
Thank you very much for your attention to this inquiry. I look forward to your response.
Best regards,
Chen
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;
}
}
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
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.
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.
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.
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.
Hello Daniel,
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 1.
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
I initialize three matrix for BFGS method :
identity = fmatrix(1, NLBFGS_vec, 1, NLBFGS_vec);
inv_hes = fmatrix(1, NLBFGS_vec, 1, NLBFGS_vec);
h_inv_plus = fmatrix(1, NLBFGS_vec, 1, NLBFGS_vec);
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.
The elastic FWI of hydrophone data (QUELLTYPB = 4) is currently not supported.
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.
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.
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
@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
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?
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.