# Creating input files for DENISE Black-Edition

Jupyter notebook for the definition of FD model and modelling/FWI/RTM parameters of DENISE Black-Edition. For a more detailed explanation of the parameters, I refer to the DENISE Black-Edition user manual

Daniel Koehn

Kiel, 24.06.2019



In :

# Import Python libaries
# ----------------------
import numpy as np                   # NumPy library
from denise_IO.denise_out import *   # "DENISE" library



## 1. Short description of modelling/FWI problem

Define name of DENISE parameter file



In :

para["filename"] = "DENISE_marm_OBC.inp"



Give a short description of your modelling/FWI problem



In :

para["descr"] = "Marmousi-II"



What kind of PHYSICS do you want to use? (2D-PSV=1; 2D-AC=2; 2D-PSV-VTI=3; 2D-PSV-TTI=4; 2D-SH=5)



In :

para["PHYSICS"] = 1



Choose DENISE operation mode (MODE): (forward_modelling_only=0; FWI=1; RTM=2)



In :

para["MODE"] = 0



## 2. Load external 2D elastic model

First define spatial model discretization:



In :

para["NX"] = 500     # number of grid points in x-direction
para["NY"] = 174     # number of grid points in y-direction
para["DH"] = 20.     # spatial grid point distance [m]





In :

# Define model basename
base_model = "model/marmousi_II_marine"

# Open vp-model and write IEEE-le binary data to vp array
# -------------------------------------------------------
f = open(base_model + ".vp")
data_type = np.dtype ('float32').newbyteorder ('<')
vp = np.fromfile (f, dtype=data_type)
f.close()

# Reshape (1 x nx*ny) vector to (ny x nx) matrix
vp = vp.reshape(para["NX"],para["NY"])
vp = np.transpose(vp)
vp = np.flipud(vp)

# Open vs-model and write IEEE-le binary data to vs array
# -------------------------------------------------------
f = open(base_model + ".vs")
data_type = np.dtype ('float32').newbyteorder ('<')
vs = np.fromfile (f, dtype=data_type)
f.close()

# Reshape (1 x nx*ny) vector to (ny x nx) matrix
vs = vs.reshape(para["NX"],para["NY"])
vs = np.transpose(vs)
vs = np.flipud(vs)

# Open rho-model and write IEEE-le binary data to rho array
# ---------------------------------------------------------
f = open(base_model + ".rho")
data_type = np.dtype ('float32').newbyteorder ('<')
rho = np.fromfile (f, dtype=data_type)
f.close()

# Reshape (1 x nx*ny) vector to (ny x nx) matrix
rho = rho.reshape(para["NX"],para["NY"])
rho = np.transpose(rho)
rho = np.flipud(rho)



Define coordinate axis



In :

x = np.arange(para["DH"], para["DH"] * (para["NX"] + 1), para["DH"])
y = np.arange(para["DH"], para["DH"] * (para["NY"] + 1), para["DH"])

# convert m -> km
x = np.divide(x,1000.0);
y = np.divide(y,1000.0);



Plot external model



In :

cmap = "magma"  # colormap

# define minimum and maximum material parameter values
vpmin = np.min(vp)
vpmax = np.max(vp)

vsmin = np.min(vs)
vsmax = np.max(vs)

rhomin = np.min(rho)
rhomax = np.max(rho)

# plot elastic model
plot_model(vp,vs,rho,x,y,cmap,vpmin,vpmax,vsmin,vsmax,rhomin,rhomax)




C:\Users\daniel_koehn\Anaconda3\lib\site-packages\matplotlib\font_manager.py:1241: UserWarning: findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans.
(prop.get_family(), self.defaultFamily[fontext]))



Write model to IEEE-le binary file



In :

# model basename
model_basename = "marmousi_II_marine"

# location of model files
para["MFILE"] = "start/" + model_basename

# writing P-wave velocity model to IEEE-le binary file
name_model = model_basename + ".vp"
f = open (name_model, mode='wb')
data_type = np.dtype ('float32').newbyteorder ('<')
vp1 = np.array(vp, dtype=data_type)
vp1 = np.rot90(vp1,3)
vp1.tofile(f)
f.close()

# writing S-wave velocity model to IEEE-le binary file
name_model = model_basename + ".vs"
f = open (name_model, mode='wb')
data_type = np.dtype ('float32').newbyteorder ('<')
vs1 = np.array(vs, dtype=data_type)
vs1 = np.rot90(vs1,3)
vs1.tofile(f)
f.close()

# writing density model to IEEE-le binary file
name_model = model_basename + ".rho"
f = open (name_model, mode='wb')
data_type = np.dtype ('float32').newbyteorder ('<')
rho1 = np.array(rho, dtype=data_type)
rho1 = np.rot90(rho1,3)
rho1.tofile(f)
f.close()



To check if the models are correctly written to the binary files, you can use the Seismic Unix function ximage



In :

print("ximage n1=" + str(para["NY"]) + " < " + model_basename + ".vp")
print("ximage n1=" + str(para["NY"]) + " < " + model_basename + ".vs")
print("ximage n1=" + str(para["NY"]) + " < " + model_basename + ".rho")




ximage n1=174 < marmousi_II_marine.vp
ximage n1=174 < marmousi_II_marine.vs
ximage n1=174 < marmousi_II_marine.rho



## 3. Define spatial FD operator

Spatial FD operator coefficients are based on Taylor series expansion and optimized according to Holberg (1987)



In :

# Order of spatial FD operator (2, 4, 6, 8, 10, 12)
para["FD_ORDER"] = 8

# Maximum relative group velocity error E
# (minimum number of grid points per shortest wavelength is defined by FD_ORDER and E)
# values:
# 0 = Taylor coefficients
# 1 = Holberg coeff.: E = 0.1 %
# 2 =                 E = 0.5 %
# 3 =                 E = 1.0 %
# 4 =                 E = 3.0 %
para["max_relative_error"] = 3



Estimate the maximum frequency in the source wavelet, which can be modelled by the given FD grid discretization and spatial FD operator



In :

# maximum modelling frequency based on grid dispersion criterion for spatial FD operator
freqmax = calc_max_freq(vp,vs,para)




Minimum Vs:  881.0  m/s
Minimum Vp:  1500.0  m/s
no of gridpoints per minimum wavelength =  2.9
maximum source wavelet frequency =  15.189655172413794  Hz



If you want to model higher frequency wave propagation, you have to decrease the spatial gridpoint distance DH by resampling the model

## 4. Parallelization by Domain Decomposition



In :

para["NPROCX"] = 5  # number of processors in x-direction
para["NPROCY"] = 3  # number of processors in y-direction



Check if the spatial domain decomposition is consistent with the spatial FD grid discretization. The following conditions have to be satisfied

NX % NPROCX = 0

NY % NPROCY = 0



In :

check_domain_decomp(para)




Domain decomposition in x-direction is correct NX % NPROCX = 0
Domain decomposition in y-direction is correct NY % NPROCY = 0



If the domain decomposition conditions are not satisfied, you have to add additional gridpoints at the bottom and right model boundary.

## 5. Time stepping

Calculate maximum time step DT according to the Courant-Friedrichs-Lewy (CFL) criterion



In :

para["DT"] = check_stability(vp,vs,para)




Maximum Vs:  2752.0  m/s
Maximum Vp:  4766.604  m/s
According to the Courant-Friedrichs-Lewy (CFL) criterion
the maximum time step is DT = 2.14e-03 s



If you want to apply a FWI, keep in mind that the FWI will change the velocity model. Therefore, the maxium seismic velocities in the model will increase and you should choose a smaller time step than the DT derived from the CFL criterion



In :

para["TIME"] = 6.0    # time of wave propagation [s]
para["DT"] = 2.0e-3   # timestep [s]



## 6. Define acquisition geometry

### a) Receiver properites and positions

Place receivers on FD modelling grid



In :

drec =  20.                                         # receiver spacing [m]
xrec1 = 800.                                        # 1st receiver position [m]
xrec2 = 8780.                                       # last receiver position [m]
xrec = np.arange(xrec1, xrec2 + para["DH"], drec)   # receiver positions in x-direction [m]

# place receivers at depth yrec [m]
depth_rec = 460.            # receiver depth [m]
yrec = depth_rec * xrec/xrec

# assemble vectors into an array
tmp = np.zeros(xrec.size, dtype=[('var1', float), ('var2', float)])
tmp['var1'] = xrec
tmp['var2'] = yrec

# write receiver positions to file
np.savetxt(basename_rec + ".dat", tmp, fmt='%4.3f %4.3f')



Define type of seismograms SEISMO:

SEISMO=0: no seismograms

SEISMO=1: particle-velocities

SEISMO=2: pressure (hydrophones)

SEISMO=3: curl and div

SEISMO=4: everything



In :

# type of seismogram
para["SEISMO"] = 1



How does DENISE read receiver positions from a file? In case of a fixed spread geometry you only need a single receiver file (READREC=1). If you want to model streamer geometry or more generally variable acquisition geometry with changing receiver positions for each shot, you have to define a separate receiver file for each shot (READREC=2)



In :



Define location and basename of receiver file, defined above, without ".dat" extension



In :



Define the seismogram properties



In :

para["NDT"] = 1   # seismogram sampling rate in timesteps (has to be set to NDT=1 if you run FWI)

# location and name of seismogram output files in SU format

# particle velocities (if SEISMO=1 or SEISMO=4)
para["SEIS_FILE_VX"] = "su/DENISE_MARMOUSI_x.su"  # filename for vx component
para["SEIS_FILE_VY"] = "su/DENISE_MARMOUSI_y.su"  # filename for vy component

# curl and div of wavefield (if SEISMO=3 or SEISMO=4)
para["SEIS_FILE_CURL"] = "su/DENISE_MARMOUSI_rot.su"  # filename for rot_z component ~ S-wave energy
para["SEIS_FILE_DIV"] = "su/DENISE_MARMOUSI_div.su"   # filename for div component ~ P-wave energy

# pressure field (hydrophones) (if SEISMO=2 or SEISMO=4)
para["SEIS_FILE_P"] = "su/DENISE_MARMOUSI_p.su"  # filename for pressure component



### b) Source properties and positions

Distribute sources on FD modelling grid an define source properties



In :

# source x-coordinates
dsrc =  80.                                         # source spacing [m]
xsrc1 = 800.                                        # 1st source position [m]
xsrc2 = 8780.                                       # last source position [m]
xsrc = np.arange(xsrc1, xsrc2 + para["DH"], dsrc)   # source positions in x-direction [m]

# place sources at depth ysrc [m]
depth_src = 40.               # source depth [m]
ysrc = depth_src * xsrc/xsrc

# number of source positions
nshot = (int)(len(ysrc))

# z-coordinate = 0 due to 2D code [m]
zsrc = 0.0 * (xsrc / xsrc)

# time delay of source wavelet [s]
td = 0.0 * (xsrc / xsrc)

# center frequency of pre-defined source wavelet [Hz]
fc = 10.0 * (xsrc / xsrc)

# you can also use the maximum frequency computed from the grid dispersion
# criterion in section 3. based on spatial discretization and FD operator
# fc = (freqmax / 2.) * (xsrc / xsrc)

# amplitude of source wavelet [m]
amp = 1.0 * (xsrc / xsrc)

# angle of rotated source [°]
angle = 0.0 * (xsrc / xsrc)

# define source type:

# 2D PSV case
# -----------
# explosive sources (QUELLTYP=1)
# point forces in x- and y-direction (QUELLTYP=2,3)

# 2D SH case
# -----------
# point force in z-direction (QUELLTYP=1)

QUELLTYP = 1
src_type = QUELLTYP * (xsrc / xsrc)

# write source positions and properties to file
basename_src = "source_OBC_VSP.dat"

# create and open source file
fp = open(basename_src, mode='w')

# write nshot to file header
fp.write(str(nshot) + "\n")

# write source properties to file
for i in range(0,nshot):

fp.write('{:4.2f}'.format(xsrc[i]) + "\t" + '{:4.2f}'.format(zsrc[i]) + "\t" +  '{:4.2f}'.format(ysrc[i]) + "\t" + '{:1.2f}'.format(td[i]) + "\t" + '{:4.2f}'.format(fc[i]) + "\t" + '{:1.2f}'.format(amp[i]) + "\t" + '{:1.2f}'.format(angle[i]) + "\t" + str(src_type[i]) + "\t" + "\n")

# close source file
fp.close()



Define location of the source file:



In :

para["SOURCE_FILE"] = "./source/" + basename_src



Do you want to excite all source positions simultaneously (RUN_MULTIPLE_SHOTS=0) or start a separate modelling run for each shot (RUN_MULTIPLE_SHOTS=1)



In :

para["RUN_MULTIPLE_SHOTS"] = 1



Define shape of the source signal (QUELLART)

• Ricker wavelet = 1
• Fuchs-Mueller wavelet = 2
• read wavelet from ASCII file = 3

• SIN^3 wavelet = 4

• Gaussian derivative wavelet = 5
\begin{equation} \rm{gd(t)=-2 \pi^2 f_c^2 (t-t_d) exp(-\pi^2 f_c^2 (t-t_d)^2)} \label{eq_s4} \end{equation}
• Bandlimited spike wavelet = 6

• Klauder wavelet = 7

\begin{equation} \rm{klau(t) = real\biggl\{sin\biggl(\frac{\pi k \tau (TS-\tau)}{\pi k \tau}\biggr)(exp(2 \pi i f_0 \tau))\biggr\}} \quad \mbox{with} \quad \rm{\tau=(t-1.5/FC\_SPIKE\_1-t_d)}\} \label{eq_s5} \end{equation}

with

$\rm{k=(FC\_SPIKE\_2-FC\_SPIKE\_1)/TS}$ (rate of change of frequency with time)

$\rm{f_0=(FC\_SPIKE\_2+FC\_SPIKE\_1)/2}$ (midfrequency of bandwidth)

$\rm{i^2=-1}$

In these equations, t denotes time and $f_c$ is the center frequency. $t_d$ is a time delay which can be defined for each source position in SOURCE_FILE. Note that the symmetric (zero phase) Ricker signal is always delayed by $1.0/f_c$, which means that after one period the maximum amplitude is excited at the source location.



In :

para["QUELLART"] = 6



If you read the wavelet from an ASCII file (QUELLART=3), you have to define the location of the signal file (SIGNAL_FILE)



In :

para["SIGNAL_FILE"] = "./wavelet/wavelet_marmousi"



In case of the bandlimited spike wavelet you have to define ...

• If FC_SPIKE_1 <= 0.0 a low-pass filtered spike with upper corner frequency FC_SPIKE_2 and order ORDER_SPIKE is calculated

• If FC_SPIKE_1 > 0.0 a band-pass filtered spike with lower corner frequency FC_SPIKE_1 and upper corner frequency FC_SPIKE_2 with order ORDER_SPIKE is calculated



In :

para["FC_SPIKE_1"] = -5.0    # lower corner frequency [Hz]
para["FC_SPIKE_2"] = 15.0    # upper corner frequency [Hz]

# you can also use the maximum frequency computed from the grid dispersion
# criterion in section 3. based on spatial discretization and FD operator
# para["FC_SPIKE_2"] = freqmax # upper corner frequency [Hz]

para["ORDER_SPIKE"] = 5      # order of Butterworth filter



In case of the Klauder wavelet you have to define the sweep length TS



In :

para["TS"] = 8.0  # sweep length [s]



Do you want to write the source wavelet to a SU file for each shot (WRITE_STF=1)?



In :

para["WRITE_STF"] = 1



Plot acquisition geometry relative to the subsurface model. Red stars denote the source positions and cyan triangles receiver positions



In :

cmap = "inferno"
plot_acq(vp,xrec/1000,yrec/1000,xsrc/1000,ysrc/1000,x,y,cmap,vpmin,vpmax)






## 7. Q-approximation



In :

para["L"] = 0     # number of relaxation mechanisms
para["FL"] = 40.  # relaxation frequencies [Hz]



## 8. Boundary conditions

Define boundary conditions. FREE_SURF=1 activates a free-surface boundary condition at y = 0 m. PML absorbing boundries are defined by their widht FW, damping velocity DAMPING, damping frequency FPML, degree of damping profile npower and k_max_PML



In :

# free surface boundary condition
para["FREE_SURF"] = 1      # activate free surface boundary condition

# PML boundary frame
para["FW"] = 10
para["DAMPING"] = 1500.
para["FPML"] = 10.
para["npower"] = 4.
para["k_max_PML"] = 1.



## 9. Wavefield snapshots

Output of wavefield snapshots (SNAP>0):

• particle velocities: SNAP=1
• pressure field: SNAP=2
• curl and divergence energy: SNAP=3
• both particle velocities and energy : SNAP=4


In :

para["SNAP"] = 0
para["SNAP_SHOT"] = 1     # compute and  write snapshots for shot no. SNAP_SHOT

para["TSNAP1"] = 0.002    # first snapshot [s] (TSNAP1 has to fullfill the condition TSNAP1 > DT)
para["TSNAP2"] = 3.0      # first snapshot [s]
para["TSNAPINC"] = 0.06   # snapshot increment [s]

para["IDX"] = 1           # write only every IDX spatial grid point in x-direction to snapshot file
para["IDY"] = 1           # write only every IDY spatial grid point in y-direction to snapshot file

para["SNAP_FILE"] = "./snap/waveform_forward"  # location and basename of the snapshot files



## 10. Log file name



In :

para["LOG_FILE"] = "log/Marmousi.log"  # Log file name



## FWI parameters

If you only want to run FD forward modelling runs, you can neglect the parameters below, which are fixed FWI parameters

## 11. General FWI parameters



In :

para["ITERMAX"] = 600   # maximum number of TDFWI iterations at each FWI stage defined in FWI workflow file

para["DATA_DIR"] = "su/MARMOUSI_spike/DENISE_MARMOUSI"   # location and basename of field data seismograms

para["INVMAT1"] = 1 # material parameterization for FWI (Vp,Vs,rho=1/Zp,Zs,rho=2/lam,mu,rho=3)

# x-y components = 1; y-comp = 2; x-comp = 3; p-comp = 4; x-p-comp = 5; y-p-comp = 6; x-y-p-comp = 7
para["QUELLTYPB"] = 1

# Optimization method
para["GRAD_METHOD"] = 2    # PCG = 1; LBFGS = 2

# PCG_BETA (Fletcher_Reeves=1/Polak_Ribiere=2/Hestenes_Stiefel=3/Dai_Yuan=4)
para["PCG_BETA"] = 2

# store NLBFGS update during LBFGS optimization
para["NLBFGS"] = 20

# store wavefields only every DTINV time sample for gradient computation
para["DTINV"] = 3

# FWI log file location and name
para["MISFIT_LOG_FILE"] = "Marmousi_fwi_log.dat"



## 12. FWI gradient taper functions



In :

# exponent of depth scaling for preconditioning

# Circular taper around all sources (not at receiver positions)
para["SWS_TAPER_CIRCULAR_PER_SHOT"] = 0
para["SRTSHAPE"] = 1    # SRTSHAPE: 1 = error_function; 2 = log_function

# Read taper file from external file
para["SWS_TAPER_FILE"] = 0



## 13. FWI model output



In :

# model location and basename
para["INV_MODELFILE"] = "model/modelTest"

# write inverted model after each iteration (yes=1)?
# Warning: Might require a lot of disk space
para["INV_MOD_OUT"] = 0



## 14. Bound constraints

Upper and lower limits for different model parameter classes



In :

# upper limit for vp
para["VPUPPERLIM"] = 6000.

# lower limit for vp
para["VPLOWERLIM"] = 0.

# upper limit for vs
para["VSUPPERLIM"] = 4000.

# lower limit for vs
para["VSLOWERLIM"] = 0.

# upper limit for density
para["RHOUPPERLIM"] = 3000.

# lower limit for density
para["RHOLOWERLIM"] = 1000.

# upper limit for Qs
para["QSUPPERLIM"] = 100.

# lower limit for Qs
para["QSLOWERLIM"] = 10.



## 15. Step length estimation



In :

para["EPS_SCALE"] = 0.01   # initial model update during step length estimation
para["STEPMAX"] = 6        # maximum number of attemps to find a step length during line search
para["SCALEFAC"] = 2.      # scale step during line search

# evaluate objective function only for a limited number of shots
para["TESTSHOT_START"] = 25
para["TESTSHOT_END"] = 75
para["TESTSHOT_INCR"] = 10



## 16. Trace muting



In :

# Activate trace muting (yes=1)
para["TRKILL"] = 0

# Location and name of trace mute file containing muting matrix
para["TRKILL_FILE"] = "./trace_kill/trace_kill.dat"



## 17. Time damping



In :

# Basename of picked traveltimes for each shot
# Time damping parameters are defined in the DENISE
# workflow file for each FWI stage
para["PICKS_FILE"] = "./picked_times/picks_"



## 18. Create DENISE parameter file



In :

write_denise_para(para)



## Instructions for preparing and starting a modelling/FWI run with DENISE Black-Edition

(a) Move model files to the directory DENISE-Black-Edition/par/para["MFILE"]



In :

print("mv " + model_basename + ".vp DENISE-Black-Edition/par/" + para["MFILE"] + ".vp")
print("mv " + model_basename + ".vs DENISE-Black-Edition/par/" + para["MFILE"] + ".vs")
print("mv " + model_basename + ".rho DENISE-Black-Edition/par/" + para["MFILE"] + ".rho")




mv marmousi_II_marine.vp DENISE-Black-Edition/par/start/marmousi_II_marine.vp
mv marmousi_II_marine.vs DENISE-Black-Edition/par/start/marmousi_II_marine.vs
mv marmousi_II_marine.rho DENISE-Black-Edition/par/start/marmousi_II_marine.rho



You can also copy the model files to a HPC cluster using SCP.

(b) Move source file to the directory DENISE-Black-Edition/par/para["SOURCE_FILE"]



In :

print("mv " + basename_src + " DENISE-Black-Edition/par/" + para["SOURCE_FILE"][2::])




mv source_OBC_VSP.dat DENISE-Black-Edition/par/source/source_OBC_VSP.dat



(c) Move receiver file(s) to the directory DENISE-Black-Edition/par/para["REC_FILE"]



In :

print("mv " + basename_rec + ".dat DENISE-Black-Edition/par" + para["REC_FILE"][1::] + ".dat")






(d) Move DENISE parameter file to the directory DENISE-Black-Edition/par/



In :

print("mv " + para["filename"] + " DENISE-Black-Edition/par/")




mv DENISE_marm_OBC.inp DENISE-Black-Edition/par/



(e) Within the DENISE-Black-Edition/par directory you can start the DENISE modelling run with



In :

print("mpirun -np " + str(para["NPROCX"]*para["NPROCY"]) + " ../bin/denise " + para["filename"])




mpirun -np 15 ../bin/denise DENISE_marm_OBC.inp



If you want to run a FWI, you also have to define a FWI workflow file



In :

print("mpirun -np " + str(para["NPROCX"]*para["NPROCY"]) + " ../bin/denise " + para["filename"] + " FWI_workflow.inp")




mpirun -np 15 ../bin/denise DENISE_marm_OBC.inp FWI_workflow.inp