adjtomo/seisflows

how to generate SPECFEM2D model .bin files?

rmodrak opened this issue · 8 comments

rmodrak commented

1) Choose desired mesh discretization

  • Choose length of and number of elements along x-axis by modifying nx, xmin, xmax in DATA/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 in DATA/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
rmodrak commented 
! 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
  • 👍1
rmodrak commented

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.

katigeophysics commented

Thank you Ryan! This detailed procedure will help me a lot. Best Regards!

luan-th-nguyen commented

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.

rmodrak commented

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?

luan-th-nguyen commented

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.

cesarari commented

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!

rmodrak commented