Minor changes for usability of toy example

par/in_and_out/ifos3d.inp

+#-----------------------------------------------------------------
+#                  PARAMETER FILE FOR IFOS3D
+#-----------------------------------------------------------------
+#
+# Note that z denotes the vertical direction !
+#
+#-----------------------------------------------------------------------------------
+#------------------------ MODELING PARAMETERS --------------------------------------
+#-----------------------------------------------------------------------------------
+#
+#-------------- Domain Decomposition -----------------------------
+number_of_processors_in_x-direction_(NPROCX) = 2
+number_of_processors_in_y-direction_(NPROCY) = 2
+number_of_processors_in_z-direction_(NPROCZ) = 2
+#
+#-------------------- 3-D Grid -----------------------------------
+number_of_gridpoints_in_x-direction_(NX) = 160
+number_of_gridpoints_in_y-direction_(NY) = 160 
+number_of_gridpoints_in_z-direction_(NZ) = 184
+distance_between_gridpoints(in_m)_in_x-direction_(DX) = 0.8	
+distance_between_gridpoints(in_m)_in_y-direction_(DY) = 0.8
+distance_between_gridpoints(in_m)_in_z-direction_(DZ) = 0.8
+#
+# Caution: first gridpoint is not at {0.0,0.0,0.0} but at {dx,dy,dz}!
+#
+#------------------- Order of FD Operator -----------------------
+order_of_spatial_fd_operators_(FDORDER) = 4
+# possible values are 2, 4, 6, 8, 10, 12
+fd_coefficients_(Taylor=1;Holberg=2)_(FDCOEFF) = 2
+#
+#-------------------Time Stepping -------------------------------
+time_of_wave_propagation_(in_sec)_(TIME) = 0.06
+timestep_(in_seconds)_(DT) = 5.0e-5
+#
+#--------------------Source---------------------------------------
+# Shape_of_source-signal:
+(ricker=1;fumue=2;from_SIGNAL_FILE=3;SIN**3=4;deltapulse=5)_(QUELLART) = 4
+point_source_(explosive=1;force_in_x=2;in_y=3;in_z=4;custom=5)_(QUELLTYP) = 4
+# If QUELLTYP <5 the following two lines are ignored
+force_angle_between_x_y_(in_degree)_(ALPHA) = 45.0
+force_angle_between_x_z_(in_degree)_(BETA) =  45.0
+# Plane wave (PW) excitation,if PLANE_WAVE_DEPTH>0, SRCREC is treated as 0
+depth_of_PW_excitation_(no<=0)_(in_meter)_(PLANE_WAVE_DEPTH) = 0.0direct.com/
+dip_of_PW_from_vertical_(in_degrees)_(PHI) =0.0
+duration_of_source-signal_PW_(in_seconds)_(TS) = 0.0033
+# External signal input instead of QUELLART
+SIGNAL_FILE = ?
+read_source_positions_from_SOURCE_FILE_(yes=1)_(SRCREC) = 1
+SOURCE_FILE = ./sources/sources_toy.dat  
+run_multiple_shots_defined_in_SOURCE_FILE_(yes=1)_(RUN_MULTIPLE_SHOTS) = 1
+#--------------------- Model (Input) -------------------------------------
+read_model_from_MFILE(yes=1)(READMOD) = 0
+MFILE = model/toy
+#
+#---------------------Q-approximation-----------------------------
+Number_of_relaxation_mechanisms_(L) = 0
+# (L=0: elastic, FWI only tested for L=0)
+L_Relaxation_frequencies_(FL) = 1000.0 
+Tau_value_for_entire_model_(TAU) = 0.000001
+#
+#----------------------Free Surface-------------------------------
+free_surface_(yes=1)(FREE_SURF) = 0
+#
+#--------------------Absorbing Boundary---------------------------
+# exponential damping applied
+type_of_boundary_condition_(ABS_TYPE)_(PML=1/ABS=2) = 1
+width_of_absorbing_frame_(in_grid_points)_(No<=0)_(FW) = 10
+percentage_of_amplitude_decay_at_outer_edge_(DAMPING) = 8.0
+# parameter for PML:
+Dominant_frequency_(in_Hz)_(FPML) = 200.0
+P_wave_velocity_near_the_grid_boundary(in_m/s)_(VPPML) = 6200.0
+# apply_periodic_boundary_condition_at_edges_(BOUNDARY):
+(no=0)_(left/right/front/back=1) = 0
+#
+#----------------------Snapshots----------------------------------
+# currently not in use!!
+output_of_snapshots_(SNAP)(yes>0) = 0
+# output of particle velocities: SNAP=1
+# output of pressure field: SNAP=2
+# output of curl and divergence energy: SNAP=3
+# output of both particle velocities and energy : SNAP=4
+first_snapshot_(in_sec)_(TSNAP1) = 0.01
+last_snapshot_(in_sec)_(TSNAP2) = 0.24
+increment_(in_sec)_(TSNAPINC) = 0.0075
+# Note that z denotes the vertical direction !
+increment_x-direction_(IDX) = 1
+increment_y-direction_(IDY) = 1
+increment_z-direction_(IDZ) = 1
+data-format_(SNAP_FORMAT)(ASCII(2);BINARY(3)) = 4
+basic_filename_(SNAP_FILE) = ./snap/back
+#output will look like SNAP_FILE.bin.z.000
+#if SNAP = 1,2 the following line is ignored
+(SNAP_PLANE) = 1
+#output of snapshots as energy wihout sign SNAP_PLANE=1
+#energy with sign true for x-z-plane SNAP_PLANE=2
+#energy with sign true for x-y-plane SNAP_PLANE=3
+#energy with sign true for y-z-plane SNAP_PLANE=4
+#
+#----------------------Receiver-----------------------------------
+output_of_seismograms_(SEISMO) = 1
+# SEISMO=0: no seismograms
+# SEISMO=1: particle-velocities
+# SEISMO=2: pressure (hydrophones)
+# SEISMO=3: curl and div
+# SEISMO=4: everything
+# Warning: "curl" is not really curl in 3D
+read_receiver_positions_from_file_(yes=1)_(READREC) = 0
+REC_FILE = ./receiver/receiver.dat
+reference_point_for_receiver_coordinate_system_(REFREC) = 0.0 , 0.0 , 0.0
+# if READREC=1 the following three lines are ignored 
+# Note that z denotes the vertical direction !
+position_of_first_receiver_(in_m)_(XREC1,YREC1,ZREC1) =  90.0 , 90.0 , 90.0
+position_of_last_receiver_(in_m)_(XREC2,YREC2,ZREC2)  =  90.0 , 90.0 , 90.0
+distance_between_two_adjacent_receivers_(in_gridpoints)_(NGEOPH) = 1
+# (Caution: receivers outside the grid may cause surprising results!)
+#
+#-------------------- Receiver array -------------------------------
+# parameters for horizontal plane of receivers
+number_of_planes_(no<=0)_(REC_ARRAY) = 1
+depth_of_first_(upper)_plane_(in_m)_(REC_ARRAY_DEPTH) = 24.0
+vertical_distance_between_planes_(in_m)_(REC_ARRAY_DIST) = 30.0 
+distance_between_receivers_in_x-direction_(in_gridpoints)_(DRX) = 10
+distance_between_receivers_in_y-direction_(in_gridpoints)_(DRY) = 10
+#
+#--------------------  Seismogram Output-----------------------------
+samplingrate_and_timelag_(in_timesteps!)_(NDT,NDTSHIFT) = 1 , 0
+# write every (ndt)th sample, leaving ndtshift samples at the beginning out
+# default: ndt=1 and ndtshift=0. (ndt=0 is set to ndt=1, ndt<0 are set to -ndt.)
+#
+data-format_(SEIS_FORMAT[6]) = 1
+#  0: SEG-Y (ASCII-text/native 4-byte-floats (IEEE on PC)/little endian on PC)
+#  1: SU (native 4-byte-floats (IEEE on PC)/little endian on PC)
+#  2: TEXTUAL (native ASCII)
+#  3: BINARY (IEEE-4-byte-floats on PC/little endian on PC)
+#  4: SEG-Y (ASCII-text/native 4-byte-floats (IEEE on PC)/little endian on PC)
+#  5: SEG-Y (ASCII-text/IBM-4-byte-floats on PC/big endian on PC) 
+#
+basic_filename_(SEIS_FILE) = ./su/cal_toy
+#
+#------------------------ Method --------------------------------
+#
+method_(METHOD) = 0
+# 0: only forward simulation
+# 1: conjugate_gradient_FWI
+#
+#-----------------------------------------------------------------------------------
+#------------------------ INVERSION PARAMETERS --------------------------------------
+#-----------------------------------------------------------------------------------
+#
+#-------------------------In- and Output Files--------------------------------------
+gradient_filename_(GRAD_FILE) = ./grad/toy
+model_output_filename_(MOD_OUT_FILE) = ./model/toy
+observed_data_fileneame_(SEIS_OBS_FILE) = ./su_obs/obs_toy
+external_or_internal_observed_data_(EXTOBS) = 0
+inversion_parameter_file_(INV_FILE) = ./in_and_out/workflow_toy.dat
+hessian_file_(HESS_FILE) = ./hess/toy
+#
+#-------------------------General---------------------------------------------------
+minimum/maximum_iteration_number_(>0)_(ITMIN,ITMAX) = 1 , 80
+filtering_(yes=1/no=0)_(FILT) = 1
+maximum_number_frequencies_per_iteration_(NFMAX) = 5
+number_of_timestep_per_period_used_for_inversion_(TAST) = 100
+average_model_parameter_(VP0,VS0,RHO0) = 6200.0, 3600.0, 2800.0
+# velocities  in m/s, density in kg/m³
+parameter_class_weighting_factors_for_vp/vs/rho_(WEIGHT) = 1.0, 1.0, 0.0
+# choose from 1.0 (full update) to 0.0 (no update)
+#
+#------------------------Steplength estimation----------------------------------------
+number_of_shots_used_for_steplength_estimation_(NSHOTS_STEP) = 4
+initial_test_steplength_(TESTSTEP) = 0.02
+#
+#------------------------Gradient preconditioning-------------------------------------
+Type_of_preconditioning_(DAMPTYPE) = 2
+# 0: no damping
+# 1: Circular taper around receivers
+# 2: Cicular taper around sources and receivers
+# 3: Tapering of source and receiver planes
+#
+#------------------------Hessian ----------------------------------------
+Apply_Hessian_(yes=1/no=0)_(HESS) = 0
+Read_Hessian_from_file_(READ_HESS) = 0
+Part_of_receivers_used_for_Hessian_(REC_HESS) = 1 
+#(Each REC_HESS Receiver is used to calculate the Hessian, not yet implemented)
+Water_level_Hessian_for_vp/vs/rho_(WATER_HESS) = 0.0192, 0.0192, 1.0e-14
+# 
+#------------------------L-BFGS ----------------------------------------
+Apply_L_BFGS_(LBFGS) = 0

par/in_and_out/ifos3d_toy.inp

-# $Id: fdmpi.inp,v 1.1.1.1 2009/12/09 14:29:22 sjetschny Exp $ 
 #-----------------------------------------------------------------
-#                  PARAMETER FILE FOR FDMPI
+#                  PARAMETER FILE FOR IFOS3D
 #-----------------------------------------------------------------
 #
 # Note that z denotes the vertical direction !
@@ -136,11 +135,11 @@ data-format_(SEIS_FORMAT[6]) = 1
 #  4: SEG-Y (ASCII-text/native 4-byte-floats (IEEE on PC)/little endian on PC)
 #  5: SEG-Y (ASCII-text/IBM-4-byte-floats on PC/big endian on PC) 
 #
-basic_filename_(SEIS_FILE) = ./su/cal_toy
+basic_filename_(SEIS_FILE) = ./su_obs/obs_toy
 #
 #------------------------ Method --------------------------------
 #
-method_(METHOD) = 1
+method_(METHOD) = 0
 # 0: only forward simulation
 # 1: conjugate_gradient_FWI
 #
@@ -157,8 +156,8 @@ inversion_parameter_file_(INV_FILE) = ./in_and_out/workflow_toy.dat
 hessian_file_(HESS_FILE) = ./hess/toy
 #
 #-------------------------General---------------------------------------------------
-minimum/maximum_iteration_number_(>0)_(ITMIN,ITMAX) = 1 , 80
-filtering_(yes=1/no=0)_(FILT) = 1
+minimum/maximum_iteration_number_(>0)_(ITMIN,ITMAX) = 1 , 60
+filtering_(yes=1/no=0)_(FILT) = 0
 maximum_number_frequencies_per_iteration_(NFMAX) = 5
 number_of_timestep_per_period_used_for_inversion_(TAST) = 100
 average_model_parameter_(VP0,VS0,RHO0) = 6200.0, 3600.0, 2800.0
@@ -173,9 +172,9 @@ initial_test_steplength_(TESTSTEP) = 0.02
 #------------------------Gradient preconditioning-------------------------------------
 Type_of_preconditioning_(DAMPTYPE) = 2
 # 0: no damping
-# 1: Gaussian taper around receivers
-# 2: Gaussian taper around sources and receivers
-# 3: Gaussian tapering of source and receiver planes
+# 1: Circular taper around receivers
+# 2: Circular taper around sources and receivers
+# 3: Tapering of source and receiver planes
 #
 #------------------------Hessian ----------------------------------------
 Apply_Hessian_(yes=1/no=0)_(HESS) = 0

par/sources/sources.dat

+32.0		32.0		92.0		0.005		300.0		1.0
+32.0		52.8		92.0		0.005		300.0		1.0
+32.0		74.4		92.0		0.005		300.0		1.0
+32.0		96.0		92.0		0.005		300.0		1.0
+64.0		32.0		92.0		0.005		300.0		1.0
+64.0		52.8		92.0		0.005		300.0		1.0
+64.0		74.4		92.0		0.005		300.0		1.0
+64.0		96.0		92.0		0.005		300.0		1.0
+96.0		32.0		92.0		0.005		300.0		1.0
+96.0		52.8		92.0		0.005		300.0		1.0
+96.0		74.4		92.0		0.005		300.0		1.0
+96.0		96.0		92.0		0.005		300.0		1.0
+
+
+
+%
+% Definition of (distributed) source positions: 
+% position, time-delay, centre frequency, and amplitude
+%
+% (comment line is indicated by # as first character)
+%
+% Parameters for each source node (one source node per line):
+%
+%  NSRC (one separateline)
+% 	XSRC		YSRC		ZSRC		TD		FC		AMP
+%
+% Symbols: 
+% XSRC= x-coordinate of source point [meter]
+% YSRC= y-coordinate of source point [meter] 
+% ZSRC= z-coordinate of source point [meter] (vertical)
+% TD= excitation time (time-delay) for source node [s]
+% FC= centre frequency of source signal [Hz]
+% AMP= maximum amplitude of source signal
+

src/hh_toy.c

@@ -19,7 +19,7 @@ float ***  taus, float ***  taup, float *  eta){
 	float Vp, Vs, Rho;
 	float *pts=NULL, ts=0.0, tp=0.0, sumu, sumpi;
 	int i, j, k, l, ii, jj, kk;
-	int kasten=0;
+	int kasten=1;
 
 	/* parameters for layer 1 */
 	const float vp1=6200.0, vs1=3600.0, rho1=2800.0;