Skip to content

Commit

Permalink
Merge pull request SPECFEM#1675 from danielpeter/devel
Browse files Browse the repository at this point in the history
adds possibility to read external source time function files in binary format
  • Loading branch information
danielpeter committed Jan 24, 2024
2 parents e9989d4 + 5cb9c6f commit 413bbc8
Show file tree
Hide file tree
Showing 3 changed files with 285 additions and 2 deletions.
2 changes: 1 addition & 1 deletion .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ env:
# example with openmp
- TESTID=31 TESTCOV=0 TESTDIR=24 TESTFLAGS="--enable-vectorization --enable-openmp" CUDA=false

# regular element mesh
# regular element mesh w/ external source time function (STF)
- TESTID=32 TESTCOV=0 TESTDIR=25 TESTFLAGS="--enable-vectorization" CUDA=false

# adjoint simulation -- offline, checked on github
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,233 @@
#!/usr/bin/env python
#
# creates an external source time function file
#
# example:
# ./create_external_source_time_function.py OUTPUT_FILES/plot_source_time_function.txt
#
# requires access to: DATA/Par_file
#
import sys
import os
import array

import numpy as np


def write_ascii_file(filename,data):
"""
writes data array to ascii file
"""
with open(filename, mode='w') as f:
# user output
print(" filename: ",filename)
print(" number of points in array = ",len(data))

# comment lines
f.write("# Source time function\n")
f.write("# number of time steps: {}\n".format(len(data)))
# data values
f.write("# STF values\n")
for i in range(len(data)):
val = data[i]
f.write("{}\n".format(val))

print(" file written")
print("")


def write_binary_file(filename,data):
"""
writes data array to binary file
"""
# float 'f' by default (otherwise 'd' for double precision)
custom_type = 'f'

with open(filename, mode='wb') as f:
# gets array length in bytes
# marker
binlength = array.array('i')
num_points = data.size
if custom_type == 'f':
# float (4 bytes) for single precision
binlength.fromlist([num_points * 4])
else:
# double precision
binlength.fromlist([num_points * 8])

# user output
print(" filename: ",filename)
print(" array length = ",binlength," Bytes")
print(" number of points in array = ",num_points)

# writes array data
binvalues = array.array(custom_type)

data = np.reshape(data, (num_points), order='F') # fortran-like index ordering
#print("debug: data ",data.tolist())

binvalues.fromlist(data.tolist()) #.tolist())

# fortran binary file: file starts with array-size again
binlength.tofile(f)
# data
binvalues.tofile(f)
# fortran binary file: file ends with array-size again
binlength.tofile(f)

print(" file written")
print("")


def create_external_source_time_function(input_file,use_ascii=True):
"""
creates an external STF file in DATA/
"""

print("")
print("creating external source time function file: DATA/stf.dat")
print(" input file : ",input_file)
print("")

if len(input_file) == 0: usage()

# checks if file exists
if not os.path.isfile(input_file):
print("Please check if file exists: ",input_file)
sys.exit(1)

# reads in input file
with open(input_file, 'r') as f:
content = f.readlines()

# checks number of lines
if len(content) == 0:
print("File is empty, please check input file")
sys.exit(1)

# get data entries
stf_data = np.array([])
for line in content:
# count valid data lines
if "#" in line or "!" in line:
# comment line
continue
else:
# assumes this is a data line
# cleans line from whitespace and extra characters
line = line.strip()
# check for empty line
if len(line) == 0: continue

# reads in value(s)
numbers = line.split()

if len(numbers) == 0:
# nothing to read
continue
if len(numbers) == 1:
# single data entry per line
val = float(numbers[0])
elif len(numbers) == 2:
# line contains two numbers (like plot_source_time_function.txt)
# we'll take the second
val = float(numbers[1])
else:
# line contains multiple numbers
# we'll take the last
val = float(numbers[-1])

# store data value
stf_data = np.append(stf_data,val)

print("number of STF time steps found: ",len(stf_data))
print("")

# checks with Par_file setting
filename = "DATA/Par_file"
with open(filename,'r') as f:
content_par = f.readlines()

NSTEP = None
for line in content_par:
# skip comments
if "#" in line or "!" in line: continue

# search for NSTEP
if "NSTEP" in line and "=" in line:
NSTEP = int(line.split('=')[1])

print("Par_file:")
print(" simulation NSTEP = ",NSTEP)
print("")

# checks
if NSTEP is None:
print("Couldn't find NSTEP in DATA/Par_file")
sys.exit(1)

# checks lengths
# must have at least NSTEP data entries in input file
if NSTEP > len(stf_data):
print("Error: simulation length NSTEP = {} is longer than the available STF data {}".format(NSTEP,len(stf_data)))
print(" Please use another input file...")
sys.exit(1)

# warning if data length is too long
if NSTEP < len(stf_data):
print("WARNING: simulation length NSTEP == {} is shorter than the available STF data {}".format(NSTEP,len(stf_data)))
print(" Only the first {} data entries will be written to the output STF file, all the others will be ignored".format(NSTEP))
print("")

# cuts STF data to simulation length
stf_data = stf_data[0:NSTEP]

# stats
stf_min = stf_data.min()
stf_max = stf_data.max()
print("stats:")
print(" source time function values: min/max = {} / {}".format(stf_min,stf_max))
print("")

# output STF file
if use_ascii:
filename = "DATA/stf.dat"
print("writing ascii output:")
write_ascii_file(filename,stf_data)
else:
filename ="DATA/stf.dat.bin"
print("writing binary output:")
write_binary_file(filename,stf_data)

print("written external source time function to: ",filename)
print("")
print("all done")
print("")


def usage():
print("usage: ./create_external_source_time_function.py input_file [--ascii or --binary]")
print(" with")
print(" input_file - ascii file with trace, e.g., OUTPUT_FILES/plot_source_time_function.txt")
print(" --ascii - (optional) writes out STF as ascii file (default)")
print(" --binary - (optional) writes out STF as binary file")
sys.exit(1)


if __name__ == '__main__':
# init
input_file = ""
use_ascii = True
# reads arguments
#print("\nnumber of arguments: " + str(len(sys.argv)))
if len(sys.argv) < 1: usage()

input_file = sys.argv[1]
if len(sys.argv) == 3:
if sys.argv[2] == "--ascii":
use_ascii = True
elif sys.argv[2] == "--binary":
use_ascii = False

# main routine
create_external_source_time_function(input_file,use_ascii)
52 changes: 51 additions & 1 deletion src/specfem3D/read_external_stf.f90
Original file line number Diff line number Diff line change
Expand Up @@ -44,22 +44,68 @@ subroutine read_external_source_time_function(isource,user_source_time_function,
character(len=MAX_STRING_LEN),intent(in) :: external_source_time_function_filename

! local variables below
integer :: i,ier
integer :: i,ier,str_len
real(kind=CUSTOM_REAL) :: stf_val
character(len=256) :: line

! threshold value to issue a warning
real(kind=CUSTOM_REAL),parameter :: SMALL_STF_VAL = 0.001

! debug timing
!double precision :: tstart,tCPU
!double precision, external :: wtime

! debug timing
!tstart = wtime()

! saftey check
if (NSTEP /= NSTEP_STF) then
print *,'Error: external source time function should have NSTEP_STF = ',NSTEP_STF,' equal to NSTEP = ',NSTEP
stop 'Error invalid number of NSTEP_STF'
endif

! for STF file names like "***.bin" - ending with ".bin", we will read it in as a binary file instead of as an ASCII file
! to enhance the speed for reading in many STF files in case.

! check file ending
! ".bin" for binary source time function files
str_len = len_trim(external_source_time_function_filename)
if (str_len > 4) then
if (external_source_time_function_filename(str_len-3:str_len) == ".bin") then
! binary source time function file
! format: can only have numbers, no comment lines
! also, it must have the exact NSTEP length (and values) for the simulation as no further checks will be done
!
! open binary file
open(IO_STF,file=trim(external_source_time_function_filename),form='unformatted',status='old',action='read',iostat=ier)
if (ier /= 0) then
print *,'Error could not open binary external source file: ',trim(external_source_time_function_filename)
stop 'Error opening binary external source time function file'
endif
read(IO_STF) user_source_time_function(:,isource)
close(IO_STF)

! debug timing
!tCPU = wtime() - tstart
!print *,'debug: binary read source ',isource,' - elapsed time: ',sngl(tCPU),'s'

! all done, no further checking (for speed)
return
endif
endif

! clear the array for that source
user_source_time_function(:,isource) = 0._CUSTOM_REAL

! ASCII source time function file
! format: the ASCII source time function file can use a format like
! # user comment
! ! another user comment
! stf_value1
! stf_value2
! ..
! stf_valueNSTEP
!
! opens specified file
open(IO_STF,file=trim(external_source_time_function_filename),status='old',action='read',iostat=ier)
if (ier /= 0) then
Expand Down Expand Up @@ -162,5 +208,9 @@ subroutine read_external_source_time_function(isource,user_source_time_function,
! closes external STF file
close(IO_STF)

! debug timing
!tCPU = wtime() - tstart
!print *,'debug: ascii read source ',isource,' - elapsed time: ',sngl(tCPU),'s'

end subroutine read_external_source_time_function

0 comments on commit 413bbc8

Please sign in to comment.