adds Berkeley UCB source time function feature

Author: danielpeter

src/shared/broadcast_computed_parameters.f90

@@ -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
@@ -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)

src/shared/read_parameter_file.F90

@@ -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'

src/shared/shared_par.f90

@@ -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
 
 !

src/specfem3D/comp_source_time_function.f90

@@ -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
 
@@ -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

src/specfem3D/locate_sources.f90

@@ -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
@@ -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'
@@ -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

src/specfem3D/rules.mk

@@ -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 \
@@ -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
@@ -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

src/specfem3D/setup_sources_receivers.f90

@@ -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
@@ -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')
@@ -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