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 [1]:

    
# 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 [2]:

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

Give a short description of your modelling/FWI problem



In [3]:

    
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 [4]:

    
para["PHYSICS"] = 1

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



In [5]:

    
para["MODE"] = 0

2. Load external 2D elastic model

First define spatial model discretization:



In [6]:

    
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]

Load external elastic model:



In [7]:

    
# 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 [8]:

    
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 [9]:

    
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 [10]:

    
# 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 [11]:

    
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 [12]:

    
# 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 [13]:

    
# 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 [14]:

    
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 [15]:

    
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 [16]:

    
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 [17]:

    
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 [18]:

    
# 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

# write receiver positions to file
basename_rec = 'receiver_OBC'
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 [19]:

    
# 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 [20]:

    
para["READREC"] = 1

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



In [21]:

    
para["REC_FILE"] = "./receiver/" + basename_rec

Define the seismogram properties



In [22]:

    
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 [23]:

    
# 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 [24]:

    
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 [25]:

    
para["RUN_MULTIPLE_SHOTS"] = 1

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.



In [26]:

    
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 [27]:

    
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 [28]:

    
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 [29]:

    
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 [30]:

    
para["WRITE_STF"] = 1

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



In [31]:

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

7. Q-approximation



In [32]:

    
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 [33]:

    
# 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 [34]:

    
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 [35]:

    
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 [36]:

    
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)

# 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"

12. FWI gradient taper functions



In [37]:

    
# 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

13. FWI model output



In [38]:

    
# 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 [39]:

    
# 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 [40]:

    
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 [41]:

    
# 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 [42]:

    
# 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 [43]:

    
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 [44]:

    
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 [45]:

    
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 [46]:

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









    



mv receiver_OBC.dat DENISE-Black-Edition/par/receiver/receiver_OBC.dat

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



In [47]:

    
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 [48]:

    
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 [49]:

    
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