<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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
Latest update
Step1:  1. Short description of modelling/FWI problem
Define name of DENISE parameter file
Step2:  Give a short description of your modelling/FWI problem
Step3:  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)
Step4:  Choose DENISE operation mode (MODE)
Step5:  2. Load external 2D elastic model
First define spatial model discretization
Step6:  Load external elastic model
Step7:  Define coordinate axis
Step8:  Plot external model
Step9:  Write model to IEEE-le binary file
Step10:  To check if the models are correctly written to the binary files, you can use the Seismic Unix function ximage
Step11:  3. Define spatial FD operator
Spatial FD operator coefficients are based on Taylor series expansion and optimized according to Holberg (1987)
Step12:  Estimate the maximum frequency in the source wavelet, which can be modelled by the given FD grid discretization and spatial FD operator, using the grid dispersion citerion
Step13:  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
Step14:  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
Step15:  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
Step16:  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
Step17:  6. Q-approximation
Step18:  7. 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
Step19:  8. Define acquisition geometry
a) Receiver properites and positions
Place receivers on FD modelling grid
Step20:  Check if receivers are located in computational domain and not the PMLs
Step21:  Write receiver positions to file
Step22:  Define type of seismograms SEISMO
Step23:  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)
Step24:  Define location and basename of receiver file, defined above, without ".dat" extension
Step25:  Define the seismogram properties
Step26:  b) Source properties and positions
Distribute sources on FD modelling grid an define source properties
Step27:  Check if sources are located in computational domain and not the PMLs
Step28:  Write source positions to file
Step29:  Define location of the source file
Step30:  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)
Step31:  Define shape of the source signal (QUELLART)
Ricker wavelet = 1
\begin{equation}
\rm{r(\tau)=\left(1-2\tau^2\right)\exp(-\tau^2)} \quad \mbox{with} \quad \rm{\tau=\frac{\pi(t-1.5/f_c-t_d)}{1.0/f_c}}
\label{eq_ricker}
\end{equation}
Fuchs-Mueller wavelet = 2
\begin{equation}
\rm{f_m(t)=\sin(2\pi(t-t_d)f_c)-0.5\sin(4\pi(t-t_d)f_c)} \quad \mbox{if} \quad \rm{t\in[t_d,t_d+1/fc]} \quad \mbox{else} \quad \rm{fm(t)=0}
\label{eq_fm}
\end{equation}
read wavelet from ASCII file = 3
SIN^3 wavelet = 4
\begin{equation}
\rm{s3(t)=0.75 \pi f_c \sin(\pi(t+t_d)f_c)^3}\quad \mbox{if} \quad \rm{t \in[t_d,t_d+1/fc]} \quad \mbox{else} \quad \rm{s3(t)=0}
\label{eq_s3}
\end{equation}
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.
Step32:  If you read the wavelet from an ASCII file (QUELLART=3), you have to define the location of the signal file (SIGNAL_FILE)
Step33:  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
Step34:  In case of the Klauder wavelet you have to define the sweep length TS
Step35:  Do you want to write the source wavelet to a SU file for each shot (WRITE_STF=1)?
Step36:  Plot acquisition geometry relative to the subsurface model. Red stars denote the source positions and cyan triangles receiver positions
Step37:  9. Wavefield snapshots
Output of wavefield snapshots (SNAP>0)
Step38:  10. Log file name
Step39:  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
Step40:  12. FWI gradient taper functions
Step41:  13. FWI model output
Step42:  14. Bound constraints
Upper and lower limits for different model parameter classes
Step43:  15. Step length estimation
Step44:  Check step length estimation
Step45:  16. Trace muting
Step46:  17. Time damping
Step47:  18. Create DENISE parameter file
Step48:  Define FWI workflow file
If you want to run a FWI, you also have to define a FWI workflow file ...
Define name of DENISE FWI workflow file
Step49:  Create Header for DENISE FWI workflow file
Step50:  1. FWI parameters for each inversion stage
For each inversion stage copy the cell below and define the FWI parameters. For a detailed description of the FWI workflow parameters, I refer to the user manual. "Experimental" options should be used with care and are still buggy.
Step51:  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"]
Step52:  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"]
Step53:  (c) Move receiver file(s) to the directory DENISE-Black-Edition/par/para["REC_FILE"]
Step54:  (d) Move DENISE parameter file to the directory DENISE-Black-Edition/par/
Step55:  (e) Within the DENISE-Black-Edition/par directory you can start the DENISE modelling run with
Step56:  If you want to run a FWI, you also have to define a FWI workflow file ...
(f) move the DENISE workflow file to the directory DENISE-Black-Edition/par/
Step57:  and run the FWI by typing
<ASSISTANT_TASK:>
Python Code:
# Import Python libaries 
# ----------------------
import numpy as np                   # NumPy library
from denise_IO.denise_out import *   # "DENISE" library
Explanation: 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
Latest update: 26.08.2020
End of explanation
para["filename"] = "DENISE_marm_OBC.inp"
Explanation: 1. Short description of modelling/FWI problem
Define name of DENISE parameter file
End of explanation
para["descr"] = "Marmousi-II"
Explanation: Give a short description of your modelling/FWI problem
End of explanation
para["PHYSICS"] = 1
Explanation: 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)
End of explanation
para["MODE"] = 0
Explanation: Choose DENISE operation mode (MODE): (forward_modelling_only=0; FWI=1; RTM=2)
End of explanation
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]
Explanation: 2. Load external 2D elastic model
First define spatial model discretization:
End of explanation
# 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)
Explanation: Load external elastic model:
End of explanation
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);
Explanation: Define coordinate axis
End of explanation
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)
Explanation: Plot external model
End of explanation
# model basename
model_basename = "marmousi_II_marine"
# location of model files during DENISE forward modelling run
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()
Explanation: Write model to IEEE-le binary file
End of explanation
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")
Explanation: To check if the models are correctly written to the binary files, you can use the Seismic Unix function ximage
End of explanation
# 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
Explanation: 3. Define spatial FD operator
Spatial FD operator coefficients are based on Taylor series expansion and optimized according to Holberg (1987)
End of explanation
# maximum modelling frequency based on grid dispersion criterion for spatial FD operator
freqmax = calc_max_freq(vp,vs,para)
Explanation: Estimate the maximum frequency in the source wavelet, which can be modelled by the given FD grid discretization and spatial FD operator, using the grid dispersion citerion
End of explanation
para["NPROCX"] = 5  # number of processors in x-direction
para["NPROCY"] = 3  # number of processors in y-direction
Explanation: 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
End of explanation
check_domain_decomp(para)
Explanation: 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
End of explanation
para["DT"] = check_stability(vp,vs,para)
Explanation: 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
End of explanation
para["TIME"] = 6.0    # time of wave propagation [s]
para["DT"] = 2.0e-3   # timestep [s]
Explanation: 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
End of explanation
para["L"] = 0     # number of relaxation mechanisms
para["FL"] = 40.  # relaxation frequencies [Hz]    
Explanation: 6. Q-approximation
End of explanation
# 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.
Explanation: 7. 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
End of explanation
# receiver x-coordinates
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
Explanation: 8. Define acquisition geometry
a) Receiver properites and positions
Place receivers on FD modelling grid
End of explanation
check_src_rec_pml(xrec,yrec,para,1)
Explanation: Check if receivers are located in computational domain and not the PMLs
End of explanation
# write receiver positions to file
basename_rec = 'receiver_OBC'
np.savetxt(basename_rec + ".dat", tmp, fmt='%4.3f %4.3f')
Explanation: Write receiver positions to file
End of explanation
# type of seismogram
para["SEISMO"] = 1
Explanation: 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
End of explanation
para["READREC"] = 1
Explanation: 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)
End of explanation
para["REC_FILE"] = "./receiver/" + basename_rec
Explanation: Define location and basename of receiver file, defined above, without ".dat" extension
End of explanation
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
Explanation: Define the seismogram properties
End of explanation
# 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)
Explanation: b) Source properties and positions
Distribute sources on FD modelling grid an define source properties
End of explanation
check_src_rec_pml(xsrc,ysrc,para,2)
Explanation: Check if sources are located in computational domain and not the PMLs
End of explanation
# 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()
Explanation: Write source positions to file
End of explanation
para["SOURCE_FILE"] = "./source/" + basename_src
Explanation: Define location of the source file:
End of explanation
para["RUN_MULTIPLE_SHOTS"] = 1
Explanation: 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)
End of explanation
para["QUELLART"] = 6
Explanation: Define shape of the source signal (QUELLART)
Ricker wavelet = 1
\begin{equation}
\rm{r(\tau)=\left(1-2\tau^2\right)\exp(-\tau^2)} \quad \mbox{with} \quad \rm{\tau=\frac{\pi(t-1.5/f_c-t_d)}{1.0/f_c}}
\label{eq_ricker}
\end{equation}
Fuchs-Mueller wavelet = 2
\begin{equation}
\rm{f_m(t)=\sin(2\pi(t-t_d)f_c)-0.5\sin(4\pi(t-t_d)f_c)} \quad \mbox{if} \quad \rm{t\in[t_d,t_d+1/fc]} \quad \mbox{else} \quad \rm{fm(t)=0}
\label{eq_fm}
\end{equation}
read wavelet from ASCII file = 3
SIN^3 wavelet = 4
\begin{equation}
\rm{s3(t)=0.75 \pi f_c \sin(\pi(t+t_d)f_c)^3}\quad \mbox{if} \quad \rm{t \in[t_d,t_d+1/fc]} \quad \mbox{else} \quad \rm{s3(t)=0}
\label{eq_s3}
\end{equation}
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.
End of explanation
para["SIGNAL_FILE"] = "./wavelet/wavelet_marmousi"
Explanation: If you read the wavelet from an ASCII file (QUELLART=3), you have to define the location of the signal file (SIGNAL_FILE)
End of explanation
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
Explanation: 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
End of explanation
para["TS"] = 8.0  # sweep length [s]
Explanation: In case of the Klauder wavelet you have to define the sweep length TS
End of explanation
para["WRITE_STF"] = 1
Explanation: Do you want to write the source wavelet to a SU file for each shot (WRITE_STF=1)?
End of explanation
cmap = "inferno"
plot_acq(vp,xrec/1000,yrec/1000,xsrc/1000,ysrc/1000,x,y,cmap,vpmin,vpmax)
Explanation: Plot acquisition geometry relative to the subsurface model. Red stars denote the source positions and cyan triangles receiver positions
End of explanation
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
Explanation: 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
End of explanation
para["LOG_FILE"] = "log/Marmousi.log"  # Log file name
Explanation: 10. Log file name
End of explanation
para["ITERMAX"] = 600   # maximum number of TDFWI iterations at each FWI stage defined in FWI workflow file
para["JACOBIAN"] = "jacobian/gradient_Test"              # location and basename of FWI gradients 
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)
# Currently, only the Vp-Vs-rho parametrization (INVMAT1=1) can be used
para["GRAD_FORM"] = 1 # gradient formulation (time integration of adjoint sources = 1, no time integration = 2)
# Adjoint source type
# 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"
Explanation: 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
End of explanation
# gradient taper geometry
para["GRADT1"] = 21
para["GRADT2"] = 25
para["GRADT3"] = 490
para["GRADT4"] = 500
para["TAPERLENGTH"] = (int)(para["GRADT2"]-para["GRADT1"])
# apply vertical taper (SWS_TAPER_GRAD_VERT=1)
para["SWS_TAPER_GRAD_VERT"] = 0
# apply horizontal taper (SWS_TAPER_GRAD_HOR=1)
para["SWS_TAPER_GRAD_HOR"] = 1
# exponent of depth scaling for preconditioning
para["EXP_TAPER_GRAD_HOR"] = 2.0
# Circular taper around all sources (not at receiver positions)
para["SWS_TAPER_GRAD_SOURCES"] = 0
para["SWS_TAPER_CIRCULAR_PER_SHOT"] = 0 
para["SRTSHAPE"] = 1    # SRTSHAPE: 1 = error_function; 2 = log_function
para["SRTRADIUS"] = 5.  # --> minimum for SRTRADIUS is 5x5 gridpoints
# Read taper file from external file
para["SWS_TAPER_FILE"] = 0
# Location and basename of taper files
para["TFILE"] = "taper/taper"
Explanation: 12. FWI gradient taper functions
End of explanation
# 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
Explanation: 13. FWI model output
End of explanation
# 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.
Explanation: 14. Bound constraints
Upper and lower limits for different model parameter classes
End of explanation
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
Explanation: 15. Step length estimation
End of explanation
check_steplength(nshot,para)
Explanation: Check step length estimation
End of explanation
# 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"
Explanation: 16. Trace muting
End of explanation
# 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_"
Explanation: 17. Time damping
End of explanation
write_denise_para(para)
Explanation: 18. Create DENISE parameter file
End of explanation
para["filename_workflow"] = "FWI_workflow_marmousi.inp"
Explanation: Define FWI workflow file
If you want to run a FWI, you also have to define a FWI workflow file ...
Define name of DENISE FWI workflow file
End of explanation
write_denise_workflow_header(para)
Explanation: Create Header for DENISE FWI workflow file
End of explanation
# Define FWI parameters for stage 1 ...
# Termination criterion
para["PRO"] = 0.01
# Frequency filtering
# TIME_FILT = 0 (apply no frequency filter to field data and source wavelet)
# TIME_FILT = 1 (apply low-pass filter to field data and source wavelet)
# TIME_FILT = 2 (apply band-pass filter to field data and source wavelet)
para["TIME_FILT"] = 1
# Low- (FC_LOW) and high-pass (FC_HIGH) corner frequencies of Butterwortfilter
# of order ORDER
para["FC_LOW"] = 0.0
para["FC_HIGH"] = 2.0
para["ORDER"] = 6
# Time windowing
para["TIME_WIN"] = 0
para["GAMMA"] = 20.0
para["TWIN-"] = 0.0
para["TWIN+"] = 0.0
# Starting FWI of parameter class Vp, Vs, rho, Qs from iteration number 
# INV_VP_ITER, INV_VS_ITER, INV_RHO_ITER, INV_QS_ITER
para["INV_VP_ITER"] = 0
para["INV_VS_ITER"] = 0
para["INV_RHO_ITER"] = 0
para["INV_QS_ITER"] = 0
# Apply spatial Gaussian filter to gradients 
# SPATFILTER = 0 (apply no filter)
# SPATFILTER = 4 (Anisotropic Gaussian filter with half-width adapted to the local wavelength)
para["SPATFILTER"] = 0
# If Gaussian filter (SPATFILTER=4), define the fraction of the local wavelength in ...
# x-direction WD_DAMP and y-direction WD_DAMP1 used to define the half-width of the 
# Gaussian filter
para["WD_DAMP"] = 0.5
para["WD_DAMP1"] = 0.5
# Preconditioning of the gradient directions
# EPRECOND = 0 - no preconditioning
# EPRECOND = 1 - approximation of the Pseudo-Hessian (Shin et al. 2001)
# EPRECOND = 3 - Hessian approximation according to Plessix & Mulder (2004)
para["EPRECOND"] = 3
# Define objective function
# LNORM = 2 - L2 norm
# LNORM = 5 - global correlation norm (Choi & Alkhalifah 2012)
# LNORM = 6 - envelope objective functions after Chi, Dong and Liu (2014) - EXPERIMENTAL
# LNORM = 7 - NIM objective function after Chauris et al. (2012) and Tejero et al. (2015) - EXPERIMENTAL
para["LNORM"] = 2
# Activate Random Objective Waveform Inversion (ROWI, Pan & Gao 2020)
# ROWI = 0 - no ROWI
# ROWI = 1 - 50% GCN l2 norm / 50% AGC l2 norm (AC, PSV, SH modules only)
para["ROWI"] = 0
# Source wavelet inversion
# STF = 0 - no source wavelet inversion
# STF = 1 - estimate source wavelet by stabilized Wiener Deconvolution 
para["STF"] = 0
# If OFFSETC_STF > 0, limit source wavelet inversion to maximum offsets OFFSETC_STF
para["OFFSETC_STF"] = -4.0
# Source wavelet inversion stabilization term to avoid division by zero in Wiener Deco
para["EPS_STF"] = 1e-1
# Apply Offset mute to field and modelled seismograms
# OFFSET_MUTE = 0 - no offset mute
# OFFSET_MUTE = 1 - mute far-offset data for offset >= OFFSETC
# OFFSET_MUTE = 1 - mute near-offset data for offset <= OFFSETC
para["OFFSET_MUTE"] = 0
para["OFFSETC"] = 10
# Scale density and Qs updates during multiparameter FWI by factors 
# SCALERHO and SCALEQS, respectively
para["SCALERHO"] = 0.5
para["SCALEQS"] = 1.0
# If LNORM = 6, define type of envelope objective function (EXPERIMENTAL)
# ENV = 1 - L2 envelope objective function
# ENV = 2 - Log L2 envelope objective function
para["ENV"] = 1
# Integrate synthetic and modelled data NORDER times (EXPERIMENTAL)
para["N_ORDER"] = 0
# Write parameters to DENISE workflow file
write_denise_workflow(para)
# Define FWI parameters for stage 2 ...
# Termination criterion
para["PRO"] = 0.01
# Frequency filtering
# TIME_FILT = 0 (apply no frequency filter to field data and source wavelet)
# TIME_FILT = 1 (apply low-pass filter to field data and source wavelet)
# TIME_FILT = 2 (apply band-pass filter to field data and source wavelet)
para["TIME_FILT"] = 1
# Low- (FC_LOW) and high-pass (FC_HIGH) corner frequencies of Butterwortfilter
# of order ORDER
para["FC_LOW"] = 0.0
para["FC_HIGH"] = 5.0
para["ORDER"] = 6
# Time windowing
para["TIME_WIN"] = 0
para["GAMMA"] = 20.0
para["TWIN-"] = 0.0
para["TWIN+"] = 0.0
# Starting FWI of parameter class Vp, Vs, rho, Qs from iteration number 
# INV_VP_ITER, INV_VS_ITER, INV_RHO_ITER, INV_QS_ITER
para["INV_VP_ITER"] = 0
para["INV_VS_ITER"] = 0
para["INV_RHO_ITER"] = 0
para["INV_QS_ITER"] = 0
# Apply spatial Gaussian filter to gradients 
# SPATFILTER = 0 (apply no filter)
# SPATFILTER = 4 (Anisotropic Gaussian filter with half-width adapted to the local wavelength)
para["SPATFILTER"] = 0
# If Gaussian filter (SPATFILTER=4), define the fraction of the local wavelength in ...
# x-direction WD_DAMP and y-direction WD_DAMP1 used to define the half-width of the 
# Gaussian filter
para["WD_DAMP"] = 0.5
para["WD_DAMP1"] = 0.5
# Preconditioning of the gradient directions
# EPRECOND = 0 - no preconditioning
# EPRECOND = 1 - approximation of the Pseudo-Hessian (Shin et al. 2001)
# EPRECOND = 3 - Hessian approximation according to Plessix & Mulder (2004)
para["EPRECOND"] = 3
# Define objective function
# LNORM = 2 - L2 norm
# LNORM = 5 - global correlation norm (Choi & Alkhalifah 2012)
# LNORM = 6 - envelope objective functions after Chi, Dong and Liu (2014) - EXPERIMENTAL
# LNORM = 7 - NIM objective function after Chauris et al. (2012) and Tejero et al. (2015) - EXPERIMENTAL
para["LNORM"] = 2
# Activate Random Objective Waveform Inversion (ROWI, Pan & Gao 2020)
# ROWI = 0 - no ROWI
# ROWI = 1 - 50% GCN l2 norm / 50% AGC l2 norm (AC, PSV, SH modules only)
para["ROWI"] = 0
# Source wavelet inversion
# STF = 0 - no source wavelet inversion
# STF = 1 - estimate source wavelet by stabilized Wiener Deconvolution 
para["STF"] = 0
# If OFFSETC_STF > 0, limit source wavelet inversion to maximum offsets OFFSETC_STF
para["OFFSETC_STF"] = -4.0
# Source wavelet inversion stabilization term to avoid division by zero in Wiener Deco
para["EPS_STF"] = 1e-1
# Apply Offset mute to field and modelled seismograms
# OFFSET_MUTE = 0 - no offset mute
# OFFSET_MUTE = 1 - mute far-offset data for offset >= OFFSETC
# OFFSET_MUTE = 1 - mute near-offset data for offset <= OFFSETC
para["OFFSET_MUTE"] = 0
para["OFFSETC"] = 10
# Scale density and Qs updates during multiparameter FWI by factors 
# SCALERHO and SCALEQS, respectively
para["SCALERHO"] = 0.5
para["SCALEQS"] = 1.0
# If LNORM = 6, define type of envelope objective function (EXPERIMENTAL)
# ENV = 1 - L2 envelope objective function
# ENV = 2 - Log L2 envelope objective function
para["ENV"] = 1
# Integrate synthetic and modelled data NORDER times (EXPERIMENTAL)
para["N_ORDER"] = 0
# Write parameters to DENISE workflow file
write_denise_workflow(para)
# Define FWI parameters for stage 3 ...
# Termination criterion
para["PRO"] = 0.01
# Frequency filtering
# TIME_FILT = 0 (apply no frequency filter to field data and source wavelet)
# TIME_FILT = 1 (apply low-pass filter to field data and source wavelet)
# TIME_FILT = 2 (apply band-pass filter to field data and source wavelet)
para["TIME_FILT"] = 1
# Low- (FC_LOW) and high-pass (FC_HIGH) corner frequencies of Butterwortfilter
# of order ORDER
para["FC_LOW"] = 0.0
para["FC_HIGH"] = 10.0
para["ORDER"] = 6
# Time windowing
para["TIME_WIN"] = 0
para["GAMMA"] = 20.0
para["TWIN-"] = 0.0
para["TWIN+"] = 0.0
# Starting FWI of parameter class Vp, Vs, rho, Qs from iteration number 
# INV_VP_ITER, INV_VS_ITER, INV_RHO_ITER, INV_QS_ITER
para["INV_VP_ITER"] = 0
para["INV_VS_ITER"] = 0
para["INV_RHO_ITER"] = 0
para["INV_QS_ITER"] = 0
# Apply spatial Gaussian filter to gradients 
# SPATFILTER = 0 (apply no filter)
# SPATFILTER = 4 (Anisotropic Gaussian filter with half-width adapted to the local wavelength)
para["SPATFILTER"] = 0
# If Gaussian filter (SPATFILTER=4), define the fraction of the local wavelength in ...
# x-direction WD_DAMP and y-direction WD_DAMP1 used to define the half-width of the 
# Gaussian filter
para["WD_DAMP"] = 0.5
para["WD_DAMP1"] = 0.5
# Preconditioning of the gradient directions
# EPRECOND = 0 - no preconditioning
# EPRECOND = 1 - approximation of the Pseudo-Hessian (Shin et al. 2001)
# EPRECOND = 3 - Hessian approximation according to Plessix & Mulder (2004)
para["EPRECOND"] = 3
# Define objective function
# LNORM = 2 - L2 norm
# LNORM = 5 - global correlation norm (Choi & Alkhalifah 2012)
# LNORM = 6 - envelope objective functions after Chi, Dong and Liu (2014) - EXPERIMENTAL
# LNORM = 7 - NIM objective function after Chauris et al. (2012) and Tejero et al. (2015) - EXPERIMENTAL
para["LNORM"] = 2
# Activate Random Objective Waveform Inversion (ROWI, Pan & Gao 2020)
# ROWI = 0 - no ROWI
# ROWI = 1 - 50% GCN l2 norm / 50% AGC l2 norm (AC, PSV, SH modules only)
para["ROWI"] = 0
# Source wavelet inversion
# STF = 0 - no source wavelet inversion
# STF = 1 - estimate source wavelet by stabilized Wiener Deconvolution 
para["STF"] = 0
# If OFFSETC_STF > 0, limit source wavelet inversion to maximum offsets OFFSETC_STF
para["OFFSETC_STF"] = -4.0
# Source wavelet inversion stabilization term to avoid division by zero in Wiener Deco
para["EPS_STF"] = 1e-1
# Apply Offset mute to field and modelled seismograms
# OFFSET_MUTE = 0 - no offset mute
# OFFSET_MUTE = 1 - mute far-offset data for offset >= OFFSETC
# OFFSET_MUTE = 1 - mute near-offset data for offset <= OFFSETC
para["OFFSET_MUTE"] = 0
para["OFFSETC"] = 10
# Scale density and Qs updates during multiparameter FWI by factors 
# SCALERHO and SCALEQS, respectively
para["SCALERHO"] = 0.5
para["SCALEQS"] = 1.0
# If LNORM = 6, define type of envelope objective function (EXPERIMENTAL)
# ENV = 1 - L2 envelope objective function
# ENV = 2 - Log L2 envelope objective function
para["ENV"] = 1
# Integrate synthetic and modelled data NORDER times (EXPERIMENTAL)
para["N_ORDER"] = 0
# Write parameters to DENISE workflow file
write_denise_workflow(para)
# Define FWI parameters for stage 4 ...
# Termination criterion
para["PRO"] = 0.01
# Frequency filtering
# TIME_FILT = 0 (apply no frequency filter to field data and source wavelet)
# TIME_FILT = 1 (apply low-pass filter to field data and source wavelet)
# TIME_FILT = 2 (apply band-pass filter to field data and source wavelet)
para["TIME_FILT"] = 1
# Low- (FC_LOW) and high-pass (FC_HIGH) corner frequencies of Butterwortfilter
# of order ORDER
para["FC_LOW"] = 0.0
para["FC_HIGH"] = 20.0
para["ORDER"] = 6
# Time windowing
para["TIME_WIN"] = 0
para["GAMMA"] = 20.0
para["TWIN-"] = 0.0
para["TWIN+"] = 0.0
# Starting FWI of parameter class Vp, Vs, rho, Qs from iteration number 
# INV_VP_ITER, INV_VS_ITER, INV_RHO_ITER, INV_QS_ITER
para["INV_VP_ITER"] = 0
para["INV_VS_ITER"] = 0
para["INV_RHO_ITER"] = 0
para["INV_QS_ITER"] = 0
# Apply spatial Gaussian filter to gradients 
# SPATFILTER = 0 (apply no filter)
# SPATFILTER = 4 (Anisotropic Gaussian filter with half-width adapted to the local wavelength)
para["SPATFILTER"] = 0
# If Gaussian filter (SPATFILTER=4), define the fraction of the local wavelength in ...
# x-direction WD_DAMP and y-direction WD_DAMP1 used to define the half-width of the 
# Gaussian filter
para["WD_DAMP"] = 0.5
para["WD_DAMP1"] = 0.5
# Preconditioning of the gradient directions
# EPRECOND = 0 - no preconditioning
# EPRECOND = 1 - approximation of the Pseudo-Hessian (Shin et al. 2001)
# EPRECOND = 3 - Hessian approximation according to Plessix & Mulder (2004)
para["EPRECOND"] = 3
# Define objective function
# LNORM = 2 - L2 norm
# LNORM = 5 - global correlation norm (Choi & Alkhalifah 2012)
# LNORM = 6 - envelope objective functions after Chi, Dong and Liu (2014) - EXPERIMENTAL
# LNORM = 7 - NIM objective function after Chauris et al. (2012) and Tejero et al. (2015) - EXPERIMENTAL
para["LNORM"] = 2
# Activate Random Objective Waveform Inversion (ROWI, Pan & Gao 2020)
# ROWI = 0 - no ROWI
# ROWI = 1 - 50% GCN l2 norm / 50% AGC l2 norm (AC, PSV, SH modules only)
para["ROWI"] = 0
# Source wavelet inversion
# STF = 0 - no source wavelet inversion
# STF = 1 - estimate source wavelet by stabilized Wiener Deconvolution 
para["STF"] = 0
# If OFFSETC_STF > 0, limit source wavelet inversion to maximum offsets OFFSETC_STF
para["OFFSETC_STF"] = -4.0
# Source wavelet inversion stabilization term to avoid division by zero in Wiener Deco
para["EPS_STF"] = 1e-1
# Apply Offset mute to field and modelled seismograms
# OFFSET_MUTE = 0 - no offset mute
# OFFSET_MUTE = 1 - mute far-offset data for offset >= OFFSETC
# OFFSET_MUTE = 1 - mute near-offset data for offset <= OFFSETC
para["OFFSET_MUTE"] = 0
para["OFFSETC"] = 10
# Scale density and Qs updates during multiparameter FWI by factors 
# SCALERHO and SCALEQS, respectively
para["SCALERHO"] = 0.5
para["SCALEQS"] = 1.0
# If LNORM = 6, define type of envelope objective function (EXPERIMENTAL)
# ENV = 1 - L2 envelope objective function
# ENV = 2 - Log L2 envelope objective function
para["ENV"] = 1
# Integrate synthetic and modelled data NORDER times (EXPERIMENTAL)
para["N_ORDER"] = 0
# Write parameters to DENISE workflow file
write_denise_workflow(para)
Explanation: 1. FWI parameters for each inversion stage
For each inversion stage copy the cell below and define the FWI parameters. For a detailed description of the FWI workflow parameters, I refer to the user manual. "Experimental" options should be used with care and are still buggy.
End of explanation
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")
Explanation: 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"]
End of explanation
print("mv " + basename_src + " DENISE-Black-Edition/par/" + para["SOURCE_FILE"][2::])
Explanation: 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"]
End of explanation
print("mv " + basename_rec + ".dat DENISE-Black-Edition/par" + para["REC_FILE"][1::] + ".dat")
Explanation: (c) Move receiver file(s) to the directory DENISE-Black-Edition/par/para["REC_FILE"]
End of explanation
print("mv " + para["filename"] + " DENISE-Black-Edition/par/")
Explanation: (d) Move DENISE parameter file to the directory DENISE-Black-Edition/par/
End of explanation
print("mpirun -np " + str(para["NPROCX"]*para["NPROCY"]) + " ../bin/denise " + para["filename"])
Explanation: (e) Within the DENISE-Black-Edition/par directory you can start the DENISE modelling run with
End of explanation
print("mv " + para["filename_workflow"] + " DENISE-Black-Edition/par/")
Explanation: If you want to run a FWI, you also have to define a FWI workflow file ...
(f) move the DENISE workflow file to the directory DENISE-Black-Edition/par/
End of explanation
print("mpirun -np " + str(para["NPROCX"]*para["NPROCY"]) + " ../bin/denise " + para["filename"] + "\t" + para["filename_workflow"])
Explanation: and run the FWI by typing
End of explanation
<END_TASK>