Skip to content

GitLab

IFOS3D
IFOS3D
Commit 09394478 authored by Simone Butzer's avatar Simone Butzer
Browse files
Options

change phrases fdmpi to ifos

parent 7da2cebb
Hide whitespace changes
Inline Side-by-side
Showing with 84 additions and 1891 deletions
par/compileIFOS3D.sh
View file @ 09394478
cd ../src
make clean
make fdmpi
make ifos3d
cd ../par
par/in_and_out/fdmpi.inp deleted 100644 → 0
View file @ 7da2cebb
# $Id: fdmpi.inp,v 1.1.1.1 2009/12/09 14:29:22 sjetschny Exp $
#-----------------------------------------------------------------
# PARAMETER FILE FOR FDMPI
#-----------------------------------------------------------------
#
# Note that z denotes the vertical direction !
#
#-----------------------------------------------------------------------------------
#------------------------ MODELING PARAMETERS --------------------------------------
#-----------------------------------------------------------------------------------
#
#-------------- Domain Decomposition -----------------------------
number_of_processors_in_x-direction_(NPROCX) = 8
number_of_processors_in_y-direction_(NPROCY) = 8
number_of_processors_in_z-direction_(NPROCZ) = 8
#
#-------------------- 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) = 1
# 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
#
#------------------------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: Gaussian taper around receivers
# 2: Gaussian taper around sources and receivers
# 3: Gaussian 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/startIFOS3D.sh
View file @ 09394478
lamboot -v lamhosts
mpirun -np 8 nice -19 ../bin/fdmpi ./in_and_out/fdmpi_toy.inp | tee ./in_and_out/fdmpi.out
#lamboot -v lamhosts
mpirun -np 8 nice -19 ../bin/ifos3d ./in_and_out/ifos3d_toy.inp | tee ./in_and_out/ifos3D.out
scripts/startJURECA.sh 0 → 100644
View file @ 09394478
#!/bin/bash -l
#SBATCH --nodes=16
#SBATCH --ntasks=512
#SBATCH --ntasks-per-node=32
#SBATCH --time=12:00:00
#SBATCH --output=mpi-out.%j
#SBATCH --error=mpi_err.%j
#SBATCH --partition=batch
srun ../bin/ifos3d ./in_and_out/ifos3d_toy.inp | tee ./in_and_out/ifos3d_toy.out
src/Makefile
View file @ 09394478
......@@ -160,8 +160,8 @@ snapmerge: $(SNAPMERGE_OBJ)
# part_model: $(PARTMODEL_OBJ)
# $(CC) $(LFLAGS) $(PARTMODEL_OBJ) -o ../bin/partmodel
fdmpi: $(FDMPI_OBJ)
$(CC) $(SFLAGS) $(FDMPI_OBJ) $(LFLAGS) -o ../bin/fdmpi
ifos3d: $(FDMPI_OBJ)
$(CC) $(SFLAGS) $(FDMPI_OBJ) $(LFLAGS) -o ../bin/ifos3d
clean:
......
src/Makefile~ deleted 100644 → 0
View file @ 7da2cebb
# $Id: Makefile,v 1.1.1.1 2009/12/09 14:29:22 sjetschny Exp $
# Makefile for fdmpi
#--------------------------------------------------------
# edit here:
# source code for model generation (acoustic)
MODEL_SCR_A = hh_acoustic.c
# source code for model generation (elastic)
#MODEL_SCR = hh_bench2.c
#MODEL_SCR = hh_3DSchach.c
#MODEL_SCR = hh_FLnodes.c
#MODEL_SCR = hh_CPML.c
#MODEL_SCR = tuebbing.c
# reading input data from ASCII-file:
# MODEL_SCR = readmod.c
#MODEL_SCR = hh_toy.c
MODEL_SCR = hh_reflect.c
# Compiler (LAM: CC=hcc, CRAY T3E: CC=cc, SHARCNet: mpicc)
# ON Linux cluster running LAM
CC=hcc
LFLAGS=-lm -lmpi -lcseife
CFLAGS=-Wall -O4
SFLAGS=-L./../libcseife
IFLAGS=-I./../libcseife
#Hermit
#CC=cc
#LFLAGS=-lm -lcseife
#CFLAGS=-O3
#SFLAGS=-L./../libcseife
#IFLAGS=-I./../libcseife
#CFLAGS=-Wall -g -g3 -ggdb
# -O4
# On CRAY T3E
#CC=cc
# Juropa
#CC=mpicc
#LFLAGS=-lm
#CFLAGS=-Wall -O4 -ipo
# On HLRN system
#CC=mpcc
#LFLAGS=-lm
# ON Linux cluster ALTIX.RZ.UNI-KIEL:DE
#CC=icc
#CFLAGS=-mp -O3 -ip0q
#LFLAGS=-lmpi -lm
# ON Linux cluster ALTIX.HRZ.TU-FREIBERG.DE
#CC=icc
#CFLAGS=-mp -O3 -ip0
#LFLAGS=-lmpi -lm -i-static
# after this line, no further editing should be necessary
# --------------------------------------------------------
.c.o:
$(CC) -c $(CFLAGS) $< $(IFLAGS)
SEISMERGE_SCR = seismerge.c
SNAPMERGE_SCR = \
snapmerge.c \
merge.c \
read_par.c \
readdsk.c \
writedsk.c \
util.c
PARTMODEL_SCR = \
part_model.c \
read_par.c \
util.c
FDMPI_UTIL = \
absorb.c \
av_mat.c \
comm_ini.c \
info.c \
initproc.c \
initproc_block.c \
initsour.c \
merge.c \
mergemod.c \
note.c \
outseis.c \
output_source_signal.c \
PML_ini.c \
rd_sour.c \
read_checkpoint.c\
readdsk.c \
read_par.c \
exchange_par.c \
receiver.c \
save_checkpoint.c\
saveseis.c \
sources.c \
splitrec.c \
splitsrc.c \
timing.c \
util.c \
wavelet.c \
writedsk.c \
writemod.c \
writepar.c \
zero.c \
zero_fdmpi.c \
rwsegy.c \
FDMPI_SRC = \
fdmpi.c \
comm_ini_s.c \
checkfd_ssg.c \
CPML_coeff.c \
CPML_ini_elastic.c\
seismo_ssg.c \
matcopy.c \
surface_ssg.c \
surface_ssg_elastic.c\
update_s_ssg.c \
update_s_ssg_CPML.c \
update_s_ssg_elastic.c \
update_s_ssg_CPML_elastic.c \
update_v_ssg.c \
update_v_ssg_CPML.c \
snap_ssg.c \
exchange_v.c \
exchange_s.c \
psource.c \
readseis.c \
residual.c \
zero_wavefield.c \
gradients.c \
gradient_F.c \
conjugategrad.c \
hess_F.c \
readhess.c \
modelupdate.c \
zero_invers.c \
zero_grad.c \
precongrad.c \
hess_apply.c \
outgrad.c \
outmod.c \
steplength.c \
cpmodel.c \
readinv.c \
disc_fourier.c \
exchange_Fv.c \
model_gauss.c \
filt_seis.c \
model2_5D.c \
smooth.c \
lbfgs.c \
lbfgs_save.c \
constant_boundary.c \
$(MODEL_SCR) \
$(FDMPI_UTIL)
FDMPI_SRC_RSG = \
fdmpi_rsg.c \
checkfd_rsg.c \
seismo_rsg.c \
update_s_rsg.c \
matcopy.c \
update_v_rsg.c \
snap_rsg.c \
exchange_v_rsg.c \
exchange_s_rsg.c \
psource_rsg.c \
$(MODEL_SCR) \
$(FDMPI_UTIL)
FDMPI_ACOUSTIC_SRC = \
av_mat_acoustic.c \
fdmpi_acoustic.c \
checkfd_acoustic.c \
comm_ini_acoustic.c \
seismo_ssg_acoustic.c \
matcopy_acoustic.c \
surface_acoustic.c \
PML_ini_acoustic.c \
absorb_PML.c \
update_s_acoustic.c \
update_v_acoustic.c \
update_s_acoustic_PML.c \
update_v_acoustic_PML.c \
snap_ssg_acoustic.c \
exchange_v_acoustic.c \
exchange_s_acoustic.c \
psource_acoustic.c \
readmod_acoustic.c \
$(MODEL_SCR_A) \
$(FDMPI_UTIL)
FDMPI_OBJ = $(FDMPI_SRC:%.c=%.o)
FDMPI_OBJ_RSG = $(FDMPI_SRC_RSG:%.c=%.o)
FDMPI_ACOUSTIC_OBJ = $(FDMPI_ACOUSTIC_SRC:%.c=%.o)
SNAPMERGE_OBJ = $(SNAPMERGE_SCR:%.c=%.o)
PARTMODEL_OBJ = $(PARTMODEL_SCR:%.c=%.o)
SEISMERGE_OBJ = $(SEISMERGE_SCR:%.c=%.o)
seismerge: $(SEISMERGE_OBJ)
$(CC) $(LFLAGS) $(SEISMERGE_OBJ) -o ../bin/seismerge
snapmerge: $(SNAPMERGE_OBJ)
$(CC) $(LFLAGS) $(SNAPMERGE_OBJ) -o ../bin/snapmerge
part_model: $(PARTMODEL_OBJ)
$(CC) $(LFLAGS) $(PARTMODEL_OBJ) -o ../bin/partmodel
fdmpi: $(FDMPI_OBJ)
$(CC) $(SFLAGS) $(FDMPI_OBJ) $(LFLAGS) -o ../bin/fdmpi
fdmpi_acoustic: $(FDMPI_ACOUSTIC_OBJ)
$(CC) $(LFLAGS) $(FDMPI_ACOUSTIC_OBJ) -o ../bin/fdmpi_acoustic
clean:
find . -name "*.o" -exec rm {} \;
find . -name "*.bck" -exec rm {} \;
find . -name "core" -exec rm {} \;
src/checkfd_ssg.c
View file @ 09394478
......@@ -323,28 +323,28 @@ float *** ptaus, float *** ptaup, float *peta, float **srcpos, int nsrc, int **r
switch (SEISMO){
case 1: particle velocities only
sprintf(xfile,"%s.fdmpi.tmp",dirname(SEIS_FILE_VX));
sprintf(xfile,"%s.ifos.tmp",dirname(SEIS_FILE_VX));
if (access(xfile,W_OK|X_OK)==-1) err2(" PE %d cannot write seismograms to %s!",xfile);
sprintf(xfile,"%s.fdmpi.tmp",dirname(SEIS_FILE_VY));
sprintf(xfile,"%s.ifos.tmp",dirname(SEIS_FILE_VY));
if (access(xfile,W_OK|X_OK)==-1) err2(" PE %d cannot write seismograms to %s!",xfile);
sprintf(xfile,"%s.fdmpi.tmp",dirname(SEIS_FILE_VZ));
sprintf(xfile,"%s.ifos.tmp",dirname(SEIS_FILE_VZ));
if (access(xfile,W_OK|X_OK)==-1) err2(" PE %d cannot write seismograms to %s!",xfile);
break;
case 2 : pressure only
sprintf(xfile,"%s.fdmpi.tmp",dirname(SEIS_FILE_P));
sprintf(xfile,"%s.ifos.tmp",dirname(SEIS_FILE_P));
if (access(xfile,W_OK|X_OK)==-1) err2(" PE %d cannot write seismograms to %s!",xfile);
break;
case 4 :
sprintf(xfile,"%s.fdmpi.tmp",dirname(SEIS_FILE_VX));
sprintf(xfile,"%s.ifos.tmp",dirname(SEIS_FILE_VX));
if (access(xfile,W_OK|X_OK)==-1) err2(" PE %d cannot write seismograms to %s!",xfile);
sprintf(xfile,"%s.fdmpi.tmp",dirname(SEIS_FILE_VY));
sprintf(xfile,"%s.ifos.tmp",dirname(SEIS_FILE_VY));
if (access(xfile,W_OK|X_OK)==-1) err2(" PE %d cannot write seismograms to %s!",xfile);
sprintf(xfile,"%s.fdmpi.tmp",dirname(SEIS_FILE_VZ));
sprintf(xfile,"%s.ifos.tmp",dirname(SEIS_FILE_VZ));
if (access(xfile,W_OK|X_OK)==-1) err2(" PE %d cannot write seismograms to %s!",xfile);
case 3 : curl and div only
sprintf(xfile,"%s.fdmpi.tmp",dirname(SEIS_FILE_DIV));
sprintf(xfile,"%s.ifos.tmp",dirname(SEIS_FILE_DIV));
if (access(xfile,W_OK|X_OK)==-1) err2(" PE %d cannot write seismograms to %s!",xfile);
sprintf(xfile,"%s.fdmpi.tmp",dirname(SEIS_FILE_CURL));
sprintf(xfile,"%s.ifos.tmp",dirname(SEIS_FILE_CURL));
if (access(xfile,W_OK|X_OK)==-1) err2(" PE %d cannot write seismograms to %s!",xfile);
break;
}
......
src/fdmpi.c deleted 100644 → 0
View file @ 7da2cebb
/* $Id: fdmpi.c,v 1.1.1.1 2009/12/09 14:29:22 sjetschny Exp $ */
/*
This is program FDMPI. Viscoelastic Version.
Parallel 3-D Viscoelastic Finite Difference Seismic Modelling
T. Bohlen
PLEASE DO NOT DISTRIBUTE. PLEASE REFER OTHER PEOPLE TO THE AUTHOR.
Author
Prof. Dr. Thomas Bohlen, Karlsruhe Institute of Technology, Geophysical Institute,
Hertzstr. 16, 76227 Karlsruhe, Germany
Phone/Fax: +49 (0)721 608 4416
mailto:thomas.bohlen@kit.edu,
http://www.gpi.kit.edu
If you want to publish synthetic data calculated with this program please
give a reference to the following paper:
Bohlen, T., 2002, Parallel 3-D viscoelastic finite-recposdifference seismic modelling,
Computers @ Geopsciences, Vol. 28, No. 8, 887-889.*/
//
#include "fd.h"
#include "globvar.h"
int main(int argc, char **argv){
int ns, nt, nseismograms=0, nf1, nf2;
int lsnap, nsnap=0, lsamp=0, nlsamp=0, buffsize;
int ntr=0, ntr_loc=0, ntr_glob=0, nsrc=0, nsrc_loc=0, ishot, ishot1, nshots;
double time1=0.0, time2=0.0, time3=0.0, time4=0.0;
double * time_v_update, * time_s_update, * time_s_exchange,* time_v_exchange, * time_timestep;
int * xb, * yb, * zb, l,i,j;
float *** absorb_coeff=NULL;
float *** sxy=NULL, *** syz=NULL, *** sxz=NULL;
float *** sxx=NULL, *** syy=NULL, *** szz=NULL;
float *** rxy=NULL, *** ryz=NULL, *** rxz=NULL;
float *** rxx=NULL, *** ryy=NULL, *** rzz=NULL;
float *** vx=NULL, *** vy=NULL, *** vz=NULL;
float *** rho=NULL, *** pi=NULL, *** u=NULL;
float *** taus=NULL, *** taup=NULL, *eta=NULL;
float *** uipjp=NULL, *** ujpkp=NULL, *** uipkp=NULL, *** tausipjp=NULL, *** tausjpkp=NULL, *** tausipkp=NULL,*** rjp=NULL, *** rkp=NULL, *** rip=NULL;
float ** srcpos=NULL, **srcpos_loc=NULL, **srcpos_loc_back=NULL,** srcpos1=NULL, ** signals=NULL;
int ** recpos=NULL, ** recpos_loc=NULL, *snum_loc=NULL,*rnum_loc=NULL ;
float ** sectionvx=NULL, ** sectionvy=NULL, ** sectionvz=NULL, ** sectionp=NULL,
** sectioncurl=NULL, ** sectiondiv=NULL;
float *** bufferlef_to_rig, *** bufferrig_to_lef;
float *** buffertop_to_bot, *** bufferbot_to_top;
float *** bufferfro_to_bac, *** bufferbac_to_fro;
float *** sbufferlef_to_rig, *** sbufferrig_to_lef;
float *** sbuffertop_to_bot, *** sbufferbot_to_top;
float *** sbufferfro_to_bac, *** sbufferbac_to_fro;
float * K_x=NULL, * alpha_prime_x=NULL, * a_x=NULL, * b_x=NULL, * K_x_half=NULL, * alpha_prime_x_half=NULL, * a_x_half=NULL, * b_x_half=NULL, * K_y=NULL, * alpha_prime_y=NULL, * a_y=NULL, * b_y=NULL, * K_y_half=NULL, * alpha_prime_y_half=NULL, * a_y_half=NULL, * b_y_half=NULL, * K_z=NULL, * alpha_prime_z=NULL, * a_z=NULL, * b_z=NULL, * K_z_half=NULL, * alpha_prime_z_half=NULL, * a_z_half=NULL, * b_z_half=NULL;
float *** psi_sxx_x=NULL, *** psi_syy_y=NULL, *** psi_szz_z=NULL, *** psi_sxy_y=NULL, *** psi_sxy_x=NULL, *** psi_sxz_x=NULL, *** psi_sxz_z=NULL, *** psi_syz_y=NULL, *** psi_syz_z=NULL, *** psi_vxx=NULL, *** psi_vyy=NULL, *** psi_vzz=NULL, *** psi_vxy=NULL, *** psi_vxz=NULL, *** psi_vyx=NULL, *** psi_vyz=NULL, *** psi_vzx=NULL, *** psi_vzy=NULL;
float ** sectionread=NULL, ** sectionreadf=NULL, ** sectionvxdiff=NULL, ** sectionvydiff=NULL,** sectionvzdiff=NULL, * misfit;
float L2=0.0, L2all=0.0, L2f=0.0;
/*inversion variables*/
/*float **** fvx=NULL, **** fvy=NULL, **** fvz=NULL;*/
float *** grad1=NULL, *** grad2=NULL, *** grad3=NULL;
float *** gradprior1=NULL, *** gradprior2=NULL, *** gradprior3=NULL;
float *** gradprior4=NULL, *** gradprior5=NULL, *** gradprior6=NULL;
float *** hess1=NULL, *** hess2=NULL, *** hess3=NULL;
float **** Ffvx=NULL, **** Ffvy=NULL, **** Ffvz=NULL, **** Ffivx=NULL, **** Ffivy=NULL, **** Ffivz=NULL;
float **** Fbvx=NULL, **** Fbvy=NULL, **** Fbvz=NULL, **** Fbivx=NULL, **** Fbivy=NULL, **** Fbivz=NULL;
int iteration=0, steptest=0,cdf=0, groupnum=0, nf=0, * itpergroup,it_group=0,itmax=0, ntast=1;
int ntr_hess=0;
float * step , *finv, *beta;
float *** testrho, *** testpi, *** testu;
float buf;
/*MPI_Status status;*/
float dummy=0.0;
int hloop=0,pshot=0,pshot1=0, pshot_loc=0;
int LBFGS=0, bfgsnum=5;
float **bfgsmod1=NULL,**bfgsgrad1=NULL, *bfgsscale1=NULL; /* *bfgsscale2,*bfgsscale3;**bfgsmod2,**bfgsmod3,**bfgsgrad2, **bfgsgrad3,*/
MPI_Request *req_send, *req_rec, *sreq_send, *sreq_rec;
MPI_Status *send_statuses, *rec_statuses;
float memdyn, memmodel, memseismograms, membuffer, memtotal,memcpml=0.0,memdynf=0.0, memgrad=0.0, membfgs=0.0;
float fac1, fac2,fac3;
char *buff_addr, ext[10];
char buffer[STRING_SIZE], bufferstring[10];
/*char comp[6];*/
FILE * fpsrc=NULL;
/* Initialize MPI environment */
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&NP);
MPI_Comm_rank(MPI_COMM_WORLD,&MYID);
setvbuf(stdout, NULL, _IONBF, 0);
/* initialize clock for estimating runtime of program */
if (MYID == 0){
time1=MPI_Wtime();
clock();
}
/* print program name, version, author etc to stdout*/
if (MYID == 0) info(stdout);
FDMPIVERS=33; /* 3D isotropic elastic */
/* PE 0 is reading the parameters from the input file fdmpi.inp */
if (MYID == 0){
if (argc>1) {
strncpy(FILEINP,argv[1],STRING_SIZE);
fprintf(stderr," Input parameter filename read from command line : %s. \n\n",FILEINP);
if (strlen(FILEINP)>STRING_SIZE-1) {
fprintf(stderr,"\n FDMPI cannot handel pathes with more than %d characters.\n",STRING_SIZE-1);
fprintf(stderr," Error: FDMPI could not read input parameter file name. -> Exit. \n\n");
return -1;
}
}
else {
strcpy(FILEINP,"fdmpi.inp");
fprintf(stderr," Caution: input parameter filename set to default 'fdmpi.inp'. \n\n");
}
FP=fopen(FILEINP,"r");
read_par(FP);
}
/* PE 0 will broadcast the parameters to all others PEs */
exchange_par();
/* Print info on log-files to stdout */
if (MYID == 0) note(stdout);
/* open log-file (each PE is using different file) */
/* fp=stdout; */
/*sprintf(ext,".%i",MYID);
strcat(LOG_FILE,ext);*/
/* nodes MYIDo writes logging info to LOG_FILE or stdout */
if (MYID==0) FP=stdout; /* logging information will be written to standard output */
/* all other nodes write logging info to LOG_FILE */
if (MYID>0) {
/*if ((FP=fopen(LOG_FILE,"w"))==NULL) err(" Opening log-file failed.");
fprintf(FP," This is the log-file generated by PE %d \n\n",MYID);*/
FP=fopen("/dev/null","w");
FI=fopen("/dev/null","w"); }
fprintf(FP," This is the log-file generated by PE %d \n\n",MYID);
if (MYID==0){
FI=fopen("in_and_out/fdmpi_invers.out","w");
setvbuf(FI, NULL, _IONBF, 0);
fprintf(FI,"------------------ Inversion Parameters ----------------- \n");
}
/* domain decomposition */
initproc();
/* set some time counters */
NT=(int)ceil(TIME/DT); /* number of timesteps - replaces: NT=iround(TIME/DT); */
TIME=(NT-1)*DT; /* TIME set to true time of the last time step */
if (NDTSHIFT>NT) ns=0;