update solver

solver/Makefile

@@ -0,0 +1,51 @@
+############################################################################################
+#
+# The subroutines below are adapted from the Computational Toolkit provided in: 
+# 
+# Schuster, G. T. (2017). Seismic inversion. Society of Exploration Geophysicists.
+#
+# We kindly thank Prof. Schuster for allowing us to use these useful and efficient Fortran 
+# subroutines. Please Cite the book above if you use these subroutines.
+#
+############################################################################################
+
+
+include ./Makefile.config
+
+# level-1 modules: self-contained
+# --level-2 modules: use level-1 modules and are used by level-3 modules
+# --level-3 modules: use both level-1 and level-2 modules
+
+OBJS = $(DIR_OBJ)/global.o \
+ $(DIR_OBJ)/parser.o \
+ $(DIR_OBJ)/datatype.o \
+ $(DIR_OBJ)/mmi_mpi.o \
+ $(DIR_OBJ)/string.o \
+ $(DIR_OBJ)/su.o \
+ $(DIR_OBJ)/io.o \
+ $(DIR_OBJ)/math.o \
+ $(DIR_OBJ)/pml.o \
+ $(DIR_OBJ)/source.o \
+ $(DIR_OBJ)/fdcore.o \
+ $(DIR_OBJ)/solver.o
+
+all: dir fd2dmpi
+
+dir:
+ @-test -d $(DIR_MOD) || mkdir $(DIR_MOD)
+ @-test -d $(DIR_OBJ) || mkdir $(DIR_OBJ)
+
+fd2dmpi: $(OBJS) $(DIR_OBJ)/fd2dmpi.o
+ $(FCC) -o $@ $^ $(LDFLAGS)
+ mv $@ ../$(DIR_EXE)
+
+bin2su.x: $(OBJS) $(DIR_OBJ)/bin2su.o
+ $(FCC) -o $@ $^ $(LDFLAGS)
+ mv $@ ../$(DIR_EXE)
+
+$(DIR_OBJ)/%.o: %.f90
+ $(FCC) -o $@ $< $(F90FLAGS)
+
+clean:
+ rm -f $(DIR_EXE)/* $(DIR_OBJ)/*.o $(DIR_MOD)/*.mod;
+

solver/Makefile.config

@@ -0,0 +1,96 @@
+############################################################################################
+#
+# The subroutines below are adapted from the Computational Toolkit provided in: 
+# 
+# Schuster, G. T. (2017). Seismic inversion. Society of Exploration Geophysicists.
+#
+# We kindly thank Prof. Schuster for allowing us to use these useful and efficient Fortran 
+# subroutines. Please Cite the book above if you use these subroutines.
+#
+############################################################################################
+
+
+# Which kind of Fortran Compiler to use, gnu or ibm ?
+compiler = gnu
+
+# Compilation option
+__MPI = yes
+__OMP = no
+
+
+DIR_OBJ=.obj
+DIR_MOD=.mod
+DIR_EXE=bin
+
+
+# GNU Fortran Compiler
+ifeq (${compiler},gnu)
+ # Use MPI?
+ ifeq (${__MPI},yes)
+ FCC = mpif90
+# FCC = /usr/local/openmpi/bin/mpif90
+ else
+ FCC = f95
+ endif
+
+ # Use OpenMP?
+ ifeq (${__OMP},yes)
+ FLAGOMP = -fopenmp
+ else
+ FLAGOMP =
+ endif
+
+ # -c: compile to obj files
+ # -w: Inhibit all warning messages
+ # -O3: Level 3 Optimization
+ # -funroll-loops: Improves performance on code using iterative DO loops by unrolling them
+ # -Jdir: Specifies where to put .mod files for compiled modules (default is -J./ )
+ F90FLAGS =-O3 -w $(FLAGOMP) -c -funroll-loops -J$(DIR_MOD)
+ LDFLAGS =-O3 -w $(FLAGOMP) 
+endif
+
+
+# INTEL Fortran Compiler
+ifeq (${compiler},intel)
+ # Use MPI?
+ ifeq (${__MPI},yes)
+ FCC = mpif90
+ else
+ FCC = ifort
+ endif
+
+ # Use OpenMP?
+ ifeq (${__OMP},yes)
+ FLAGOMP = -openmp
+ else
+ FLAGOMP =
+ endif
+
+ F90FLAGS =-O3 -w $(FLAGOMP) -c -xP -fpp -fno-fnalias -fno-alias -ip -fomit-frame-pointer -ftz -assume byterecl -module $(DIR_MOD)
+ LDFLAGS =-O3 -w $(FLAGOMP) 
+endif
+
+
+# IBM Fortran Compiler
+ifeq (${compiler},ibm)
+ # Use MPI?
+ ifeq (${__MPI},yes)
+ FCC = mpixlf90
+ ifeq (${__OMP},yes)
+ FCC = mpixlf90_r
+ endif
+ else
+ FCC = f95
+ endif
+
+ # Use OpenMP?
+ ifeq (${__OMP},yes)
+ FLAGOMP = -qsmp=noauto:omp
+ else
+ FLAGOMP =
+ endif
+
+ F90FLAGS =-O3 -u -qunroll=yes $(FLAGOMP) -qarch=450d -qtune=450 -c -qmoddir=$(DIR_MOD) -I$(DIR_MOD)
+ LDFLAGS =-O3 -u -qunroll=yes $(FLAGOMP) -qarch=450d -qtune=450
+endif
+

solver/README.md

@@ -0,0 +1,10 @@
+# fd2dmpi: Acoustic Wavefield Solver
+
+#### The variable-density acoustic wave equation is solved using a time-domain staggered FD method.
+
+The Fortran subroutines are adapted from the Computational Toolkit provided in: 
+
+​ Schuster, G. T. (2017). Seismic inversion. Society of Exploration Geophysicists. https://library.seg.org/doi/book/10.1190/1.9781560803423
+
+We kindly thank Prof. Schuster for allowing us to use these useful and efficient Fortran subroutines. Please Cite the book above if you use these subroutines.
+

solver/bin2su.f90

@@ -0,0 +1,70 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! The subroutines below are adapted from the Computational Toolkit provided in: 
+! 
+! Schuster, G. T. (2017). Seismic inversion. Society of Exploration Geophysicists.
+!
+! We kindly thank Prof. Schuster for allowing us to use these useful and efficient Fortran 
+! subroutines. Please Cite the book above if you use these subroutines.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+! Program for converting a binary file into a SU format.
+!
+program bin2su
+
+use datatype
+use io
+use mmi_mpi
+use parser
+use string
+use su
+
+implicit none
+
+type(acquisition) :: coord
+character(200) :: coordfile, csg_in, data_out, input, output
+integer :: is, ig, nt_in, nt_out, iargc
+real, allocatable :: csg(:,:)
+
+if (iargc() == 0) then
+ write(*,*) 'Usage: bin2su.x coordfile csg_in data_out nt_in nt_out'
+ goto 999
+endif
+
+! Read input parameters from command line
+call readCmd('coordfile', coordfile)
+call readCmd('csg_in', csg_in)
+call readCmd('data_out', data_out)
+call readCmd('nt_in', nt_in)
+call readCmd('nt_out', nt_out)
+
+call MPI_INIT(ierr)
+call MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierr)
+call MPI_COMM_SIZE(MPI_COMM_WORLD,nsize,ierr)
+
+call readcoordfile(coordfile, coord)
+
+! Allocate memory
+allocate(csg(nt_in,coord%ngmax))
+
+! Convert data
+do is=1,coord%ns
+ write(*,*) 'is = ', is
+
+ ! Read .bin file
+ call filename(input, csg_in, is, '.bin');
+ call read_binfile(input, csg, nt_in, coord%ng(is))
+
+ ! Write .su file
+ call filename(output, data_out, is, '.su');
+ call write_sufile(output, csg(1:nt_out,1:coord%ng(is)), nt_out, coord%ng(is), 1.0, 1.0)
+enddo
+
+call MPI_Finalize(MPI_COMM_WORLD, ierr)
+
+999 continue
+
+end program bin2su
+

solver/datatype.f90

@@ -0,0 +1,39 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! The subroutines below are adapted from the Computational Toolkit provided in: 
+! 
+! Schuster, G. T. (2017). Seismic inversion. Society of Exploration Geophysicists.
+!
+! We kindly thank Prof. Schuster for allowing us to use these useful and efficient Fortran 
+! subroutines. Please Cite the book above if you use these subroutines.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+! module for defining data structures
+!
+module datatype
+
+type param
+ integer :: nx, nz, nt, free_surface, npml, store_snap, store_step
+ real :: dx, dt, f, vmin, vmax, xmin, xmax
+ character(len=200) :: coordfile, data_out, velfile, densityfile, fileformat, sourcefile, data, sourcetype
+
+end type param
+
+type acquisition
+ ! Data format: t(is,ig), xg(iz,ix), zg(iz,ix)
+
+ integer :: ns, ngmax
+ integer, allocatable :: ng(:)
+ real, allocatable :: xs(:), zs(:), xg(:,:), zg(:,:), t(:,:)
+end type acquisition
+
+type wavelet_param
+ character(len=200) :: fileformat, inputprefix, outputprefix
+ integer :: ns, ng, nt, iref
+ real :: dx, dt, ds, dg
+end type wavelet_param
+
+end module datatype
+

solver/fd2dmpi.f90

@@ -0,0 +1,43 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! The subroutines below are adapted from the Computational Toolkit provided in: 
+! 
+! Schuster, G. T. (2017). Seismic inversion. Society of Exploration Geophysicists.
+!
+! We kindly thank Prof. Schuster for allowing us to use these useful and efficient Fortran 
+! subroutines. Please Cite the book above if you use these subroutines.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+! 
+! Main program for forward and adjoint modeling.
+!
+program fd2dmpi
+
+use solver
+use parser
+
+implicit none
+
+character*256 :: parfile
+
+ ! Read the parameter file
+ call readCmd('par',parfile,'solver.config')
+ call readParFile(parfile,'jobtype',jobtype)
+
+ ! Forward modeling
+ if (jobtype(1:7) == 'forward') then
+ call forward_modeling(parfile)
+
+ ! Adjoint modeling
+ else if (jobtype(1:7) == 'adjoint') then
+ call adjoint_modeling(parfile)
+
+ ! Gradient computation
+ else if (jobtype(1:8) == 'gradient') then
+ call gradient_computing(parfile)
+ endif
+
+end program fd2dmpi
+