how to generate SPECFEM2D model .bin files?
rmodrak opened this issue · 8 comments
1) Choose desired mesh discretization
-
Choose length of and number of elements along x-axis by modifying
nx
,xmin
,xmax
inDATA/Par_file
. -
Choose length of and number of elements along z-axis by modifying
DATA/interfaces.dat
2) output SPECFEM2D mesh coordinates
-
Specify a homogenous or water-layer-over-halfspace model using
nbmodels
,nbregions
and associated parameters inDATA/Par_file
-
Choose the following settings in
DATA/Par_file
MODEL = default
SAVE_MODEL = ascii
- Run mesher and solver.
bin/xmeshfem2D
bin/xspecfem3D
3) After running the mesher and solver, there should be ascii file(s) DATA/proc******_rho_vp_vs.dat
. The first two columns of this file specify the x and z coordinates of the mesh. Interpolate your values for vp,vp,rho onto these coordinates.
4) Write interpolated values obtained from the previous step in fortran binary format. To do this you'll need a short Fortran program, such as the one in the following comment.
5) If you want, plot the newly created .bin files to make sure everything look alright (see issue #34). To use your binary files, don't forget to change back the settings in DATA/Par_file
:
MODEL = binary
SAVE_MODEL = default
! XASCII_TO_BINARY
!
! USAGE
! xascii_to_binary INPUT_FILE OUTPUT_DIR
!
program sum_kernels
implicit none
logical, parameter :: USE_RHO_VP_VS = .true.
integer, parameter :: CUSTOM_REAL = 4
integer, parameter :: MAX_STRING_LEN = 255
integer :: ipath, npath, iker, nline, i,j
real(kind=CUSTOM_REAL), allocatable, dimension(:) :: kernel1, kernel2, kernel3
real(kind=CUSTOM_REAL), allocatable, dimension(:) :: xcoord, zcoord
character(len=MAX_STRING_LEN) :: input_file, output_dir
character(len=MAX_STRING_LEN) :: arg(2), filename
integer :: ios, ier
print *, 'Running XASCII_TO_BINARY'
print *,
if (command_argument_count() /= 2) then
print *, 'USAGE: ascii_to_binary INPUT_FILE OUTPUT_DIR'
print *, ''
stop 'Please check command line arguments'
endif
! parse command line arguments
do i = 1, 2
call get_command_argument(i,arg(i), status=ier)
enddo
read(arg(1),'(a)') input_file
read(arg(2),'(a)') output_dir
call get_number_lines(input_file, nline)
print *, 'Number of lines: ', nline
print *,
allocate( xcoord(nline) )
allocate( zcoord(nline) )
allocate( kernel1(nline) )
allocate( kernel2(nline) )
allocate( kernel3(nline) )
xcoord(:) = 0.d0
zcoord(:) = 0.d0
kernel1(:) = 0.00
kernel2(:) = 0.00
kernel3(:) = 0.00
open(unit=100,file=trim(input_file),status='old',action='read')
do j=1,nline
read(100,*) xcoord(j), zcoord(j), kernel1(j), kernel2(j), kernel3(j)
enddo
close(100)
i=0
write(filename,'(a,i6.6,a)') '/proc', i, '_x.bin'
open(unit=101, file=trim(output_dir)//trim(filename),status='unknown',action='write',form='unformatted')
write(filename,'(a,i6.6,a)') '/proc', i, '_z.bin'
open(unit=102,file=trim(output_dir)//trim(filename),status='unknown',action='write',form='unformatted')
if (USE_RHO_VP_VS) then
write(filename,'(a,i6.6,a)') '/proc', i, '_rho.bin'
open(unit=103,file=trim(output_dir)//trim(filename),status='unknown',action='write',form='unformatted')
write(filename,'(a,i6.6,a)') '/proc', i, '_vp.bin'
open(unit=104,file=trim(output_dir)//trim(filename),status='unknown',action='write',form='unformatted')
write(filename,'(a,i6.6,a)') '/proc', i, '_vs.bin'
open(unit=105,file=trim(output_dir)//trim(filename),status='unknown',action='write',form='unformatted')
endif
print *, 'Writing kernels to: ', trim(output_dir)
print *,
write(101) xcoord
write(102) zcoord
write(103) kernel1
write(104) kernel2
write(105) kernel3
close(101)
close(102)
close(103)
close(104)
close(105)
print *, 'Finished writing kernels '
print *,
end program sum_kernels
! ------------------------------------------------------------------------------
subroutine get_number_lines(filename, nline)
integer, parameter :: CUSTOM_REAL = 4
integer, parameter :: MAX_STRING_LEN = 255
integer, parameter :: MAX_NUMBER_LINES = 2147483647
real(kind=CUSTOM_REAL) :: dummy1, dummy2, dummy3, dummy4, dummy5
character(len=MAX_STRING_LEN) :: filename
integer :: ios, nline
print *, trim(filename)
open(unit=100,file=trim(filename),status='old',action='read')
nline = 0
do j=1,MAX_NUMBER_LINES
read(100,*,iostat=ios) dummy1, dummy2, dummy3, dummy4, dummy5
if (ios /= 0) exit
nline=nline+1
enddo
close(100)
end subroutine get_number_lines
In my opinion, the procedure for writing SPECFEM2D models is unnecessarily complicated. Some of this complexity results from the fact that, unlike finite difference grids, spectral element meshes are not regularly spaced. But mostly this complexity results from the quick-and-dirty nature of SPECFEM2D itself!
Unfortunately, options for solvers are limited, and we are stuck using SPECFEM2D until someone writes an interface to a different open source solver.
Thank you Ryan! This detailed procedure will help me a lot. Best Regards!
I see that creating a binary model in the current version of SPECFEM2D is way easier. You just need to do until step 2 in Ryan's instructions above but with setting the parameter 'SAVE_MODEL = binary' in Par_file. The generated binary files will be placed in DATA after simulation of the default model is finished.
Hi Luan, Thanks for the update. But even with the current version of SPECFEM2D, wouldn't you still need to follow steps 3-5 if your desired model is not homogeneous?
Hi Ryan, you are right. An additional interpolation step is needed to create a model that is different from the default model defined in Par_file. Of course the default model can be homogeneous or it can contain rectangular material regions.
As far as I can remember, beside allowing to save the binary directly, the current version of SPECFEM2D also fixed the problem with reading a multi-slice model that happens in the commit #d745c542 of SPECFEM2D.
Hello. Thanks for sharing seisflow. I just downloaded and installed successfully seisflow and SPECFEM2D, but when I run
sfclean ; sfrun within ~/tests/checkers, and I get the following error:
Traceback (most recent call last):
File "/home/cesar/packages/seisflows/scripts/sfrun", line 70, in
config()
File "/home/cesar/packages/seisflows/seisflows/config.py", line 63, in config
sys.modules['seisflows_'+name] = custom_import(name)()
File "/home/cesar/packages/seisflows/seisflows/config.py", line 207, in custom_import
(args[0], args[1], args[0].upper()))
Exception:
SEISFLOWS IMPORT ERROR
The following module was not found in the SeisFlows package:
seisflows.system.chinook_sm
Please check user-supplied SYSTEM parameter.
What should I do? Thanks!