Skip to content

Commit

Permalink
adds Berkeley UCB source time function feature
Browse files Browse the repository at this point in the history
  • Loading branch information
@danielpeter
danielpeter committed Oct 10, 2024
1 parent 71decbc commit f2569f1
Show file tree
Hide file tree
Showing 8 changed files with 2,806 additions and 7 deletions.
11 changes: 10 additions & 1 deletion src/shared/broadcast_computed_parameters.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@

subroutine broadcast_computed_parameters()

use constants, only: myrank
use constants, only: myrank,STF_IS_UCB_HEAVISIDE
use shared_parameters

implicit none
Expand Down Expand Up @@ -205,6 +205,15 @@ subroutine broadcast_computed_parameters()
call bcast_all_singlei(NZ_DOUBLING_4)
call bcast_all_singlei(NZ_DOUBLING_5)

! (optional) Berkeley UCB stf
if (STF_IS_UCB_HEAVISIDE) then
call bcast_all_singledp(UCB_SOURCE_T1)
call bcast_all_singledp(UCB_SOURCE_T2)
call bcast_all_singledp(UCB_SOURCE_T3)
call bcast_all_singledp(UCB_SOURCE_T4)
call bcast_all_singledp(UCB_TAU)
endif

! (optional) scattering perturbations
call bcast_all_singlel(ADD_SCATTERING_PERTURBATIONS)
call bcast_all_singledp(SCATTERING_STRENGTH)
Expand Down
9 changes: 9 additions & 0 deletions src/shared/read_parameter_file.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -220,6 +220,15 @@ subroutine read_parameter_file()
call read_value_logical(USE_MONOCHROMATIC_CMT_SOURCE, 'USE_MONOCHROMATIC_CMT_SOURCE', ier)
if (ier /= 0) stop 'an error occurred while reading the parameter file: USE_MONOCHROMATIC_CMT_SOURCE'

! (optional) Berkeley UCB STF parameters
if (STF_IS_UCB_HEAVISIDE) then
call read_value_double_precision(UCB_SOURCE_T1,'SOURCE_T1', ier); ier = 0
call read_value_double_precision(UCB_SOURCE_T2,'SOURCE_T2', ier); ier = 0
call read_value_double_precision(UCB_SOURCE_T3,'SOURCE_T3', ier); ier = 0
call read_value_double_precision(UCB_SOURCE_T4,'SOURCE_T4', ier); ier = 0
call read_value_double_precision(UCB_TAU,'TAU', ier); ier = 0
endif

! option to save strain seismograms
call read_value_logical(SAVE_SEISMOGRAMS_STRAIN, 'SAVE_SEISMOGRAMS_STRAIN', ier)
if (ier /= 0) stop 'an error occurred while reading the parameter file: SAVE_SEISMOGRAMS_STRAIN'
Expand Down
4 changes: 4 additions & 0 deletions src/shared/shared_par.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -208,6 +208,10 @@ module shared_input_parameters
logical :: SHIFT_SIMULTANEOUS_RUNS = .false.
double precision :: FILESYSTEM_IO_BANDWIDTH = 0.d0

! UCB Source time function parameters
double precision :: UCB_SOURCE_T1 = 400.d0, UCB_SOURCE_T2 = 250.d0, UCB_SOURCE_T3 = 53.d0, UCB_SOURCE_T4 = 40.d0
double precision :: UCB_TAU = 400.d0

end module shared_input_parameters

!
Expand Down
11 changes: 9 additions & 2 deletions src/specfem3D/comp_source_time_function.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,11 @@

double precision function comp_source_time_function(t,hdur,it_index)

use constants, only: EXTERNAL_SOURCE_TIME_FUNCTION
use constants, only: EXTERNAL_SOURCE_TIME_FUNCTION, &
STF_IS_UCB_HEAVISIDE

! for berkeley source time function
use ucb_heaviside, only: comp_source_time_function_ucb_stf

implicit none

Expand All @@ -38,7 +42,10 @@ double precision function comp_source_time_function(t,hdur,it_index)
double precision, external :: comp_source_time_function_heavi
double precision, external :: comp_source_time_function_ext

if (EXTERNAL_SOURCE_TIME_FUNCTION) then
if (STF_IS_UCB_HEAVISIDE) then
! Berkeley stf
comp_source_time_function = comp_source_time_function_ucb_stf(t)
else if (EXTERNAL_SOURCE_TIME_FUNCTION) then
! external stf
comp_source_time_function = comp_source_time_function_ext(it_index)
else
Expand Down
25 changes: 22 additions & 3 deletions src/specfem3D/locate_sources.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,9 @@ subroutine locate_sources()

use specfem_par_movie, only: vtkdata_source_x,vtkdata_source_y,vtkdata_source_z

! for Berkeley stf
use shared_parameters, only: UCB_SOURCE_T1,UCB_SOURCE_T2,UCB_SOURCE_T3,UCB_SOURCE_T4,UCB_TAU

implicit none

! local parameters
Expand Down Expand Up @@ -549,8 +552,17 @@ subroutine locate_sources()
write(IMAIN,*) ' half duration in frequency: ',hdur(isource),' seconds**(-1)'
case (2)
! Heaviside
write(IMAIN,*) ' using (quasi) Heaviside source time function'
write(IMAIN,*) ' half duration: ',hdur(isource),' seconds'
if (STF_IS_UCB_HEAVISIDE) then
! Berkeley UCB stf
write(IMAIN,*) ' using Berkeley UCB (quasi) Heaviside source time function'
write(IMAIN,*) ' source T1/T2/T3/T4: ',sngl(UCB_SOURCE_T1),'/',sngl(UCB_SOURCE_T2),'/', &
sngl(UCB_SOURCE_T3),'/',sngl(UCB_SOURCE_T4)
write(IMAIN,*) ' source time-shift : ',sngl(UCB_TAU)
else
! default Heaviside
write(IMAIN,*) ' using (quasi) Heaviside source time function'
write(IMAIN,*) ' half duration: ',hdur(isource),' seconds'
endif
case (3)
! Monochromatic
write(IMAIN,*) ' using monochromatic source time function'
Expand All @@ -567,14 +579,21 @@ subroutine locate_sources()
case default
stop 'unsupported force_stf value!'
end select
else if (STF_IS_UCB_HEAVISIDE) then
! moment tensor
! Berkeley UCB stf
write(IMAIN,*) ' using Berkeley UCB (quasi) Heaviside source time function'
write(IMAIN,*) ' source T1/T2/T3/T4: ',sngl(UCB_SOURCE_T1),'/',sngl(UCB_SOURCE_T2),'/', &
sngl(UCB_SOURCE_T3),'/',sngl(UCB_SOURCE_T4)
write(IMAIN,*) ' source time-shift : ',sngl(UCB_TAU)
else if (USE_MONOCHROMATIC_CMT_SOURCE) then
! moment tensor
write(IMAIN,*) ' using monochromatic source time function'
! add message if source is monochromatic
write(IMAIN,*)
write(IMAIN,*) ' period: ',hdur(isource),' seconds'
else
! moment tensor
! default quasi Heaviside
write(IMAIN,*) ' using (quasi) Heaviside source time function'
! add message if source is a Heaviside
if (hdur(isource) <= 5.0*DT) then
Expand Down
6 changes: 6 additions & 0 deletions src/specfem3D/rules.mk
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,7 @@ specfem3D_SOLVER_OBJECTS += \
$O/setup_GLL_points.solverstatic.o \
$O/setup_sources_receivers.solverstatic.o \
$O/specfem3D_par.solverstatic_module.o \
$O/ucb_heaviside.solverstatic.o \
$O/update_displacement_LDDRK.solverstatic.o \
$O/update_displacement_Newmark.solverstatic.o \
$O/write_movie_output.solverstatic.o \
Expand Down Expand Up @@ -160,6 +161,7 @@ specfem3D_MODULES = \
$(FC_MODDIR)/siem_poisson.$(FC_MODEXT) \
$(FC_MODDIR)/siem_solver_mpi.$(FC_MODEXT) \
$(FC_MODDIR)/siem_solver_petsc.$(FC_MODEXT) \
$(FC_MODDIR)/ucb_heaviside.$(FC_MODEXT) \
$(EMPTY_MACRO)

# These files come from the shared directory
Expand Down Expand Up @@ -404,6 +406,10 @@ $O/locate_point.solverstatic.o: $O/search_kdtree.shared.o

$O/make_gravity.solver.o: $O/model_prem.shared.o $O/model_Sohl.shared.o $O/model_vpremoon.shared.o

# for Berkeley ucb stf
$O/comp_source_time_function.solverstatic.o: $O/ucb_heaviside.solverstatic.o
$O/setup_sources_receivers.solverstatic.o: $O/ucb_heaviside.solverstatic.o

$O/compute_forces_crust_mantle_Dev.solverstatic.o: $O/compute_element.solverstatic.o $O/compute_element_att_memory.solverstatic.o
$O/compute_forces_crust_mantle_noDev.solverstatic.o: $O/compute_element.solverstatic.o $O/compute_element_att_memory.solverstatic.o
$O/compute_forces_inner_core_Dev.solverstatic.o: $O/compute_element.solverstatic.o $O/compute_element_att_memory.solverstatic.o
Expand Down
16 changes: 15 additions & 1 deletion src/specfem3D/setup_sources_receivers.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -580,6 +580,10 @@ subroutine setup_sources()
use specfem_par
use specfem_par_crustmantle
use specfem_par_movie

! for berkeley ucb stf
use ucb_heaviside, only: init_ucb_heaviside

implicit none

! local parameters
Expand Down Expand Up @@ -639,11 +643,16 @@ subroutine setup_sources()
comp_dir_vect_source_Z_UP(:) = 0.d0
endif

! initialize Berkeley stf if needed (for ucb heaviside)
if (STF_IS_UCB_HEAVISIDE) then
call init_ucb_heaviside(NSTEP,DT)
endif

! sources
! moved open statement and writing of first lines into sr.vtk before the
! call to locate_sources, where further write statements to that file follow
if (myrank == 0) then
! write source and receiver VTK files for Paraview
! write source and receiver VTK files for Paraview
filename = trim(OUTPUT_FILES)//'/sr_tmp.vtk'
open(IOUT_VTK,file=trim(filename),status='unknown',iostat=ier)
if (ier /= 0 ) call exit_MPI(myrank,'Error opening temporary file sr_temp.vtk')
Expand Down Expand Up @@ -794,6 +803,11 @@ subroutine setup_stf_constants()
t0 = 0.d0
endif

! Berkeley stf
if (STF_IS_UCB_HEAVISIDE) then
t0 = 0.d0
endif

! checks if user set USER_T0 to fix simulation start time
! note: USER_T0 has to be positive
if (USER_T0 > 0.d0) then
Expand Down
Loading

0 comments on commit f2569f1

Please sign in to comment.