[CIG-SEISMO] SPECFEM3d_Cartesian , using point force with gaussian stf

[Florian Schumacher]

Hello!

I came accross the following problem with SPECFEM3d_Cartesian:


Looking at the latest code on GitHub, I noticed that the flags
USE_RICKER_TIME_FUNCTION and USE_FORCE_POINT_SOURCE can now be used
independently of each other in any combination, which allows for point
force and moment tensor sources using either, ricker or gaussian source
time functions (or error function, in case of moment tensor) .
This is great! (was strongly restricted before).


However, there seem to be some remaining old code (in devel branch)
which prevents this from working properly:

1) src/specfem3D/setup_sources_receivers.f90 , subroutine setup_sources() :
t0 is (re)defined as
  t0 = - 1.2d0 * minval(tshift_src(:) - 1.0d0/hdur(:))
in case of (USE_FORCE_POINT_SOURCE .or. USE_RICKER_TIME_FUNCTION).
This (re)definition makes sense if hdur contains dominant frequency, but
should not happen if you want to use a gaussian for a point source
(since then the array hdur should actually contain half durations).
>From what I can conclude at this point, you can fix this simply by
removing "USE_FORCE_POINT_SOURCE .or." from the if-clause, i.e. execute
this (re)definition only in case of USE_RICKER_TIME_FUNCTION == .true.

2) src/specfem3D/compute_add_sources_viscoelastic.f90:
the function comp_source_time_function_gauss() is always called with a
fixed half duration value of 5.d0*DT instead of hdur_gaussian(isource).
Changing this to hdur_gaussian(isource) would give the full
functionality that I guess is intended here (please do correct me, if
necessary).


Making these two small corrections, the correct source time functions
are created for all types of sources, having a correct time shift t0.
I have attached the two relevant files to this email, correcting the
code (see comments containing "FS FS"). However, in
compute_add_sources_viscoelastic.f90 I am not entirely sure whether all
calls to comp_source_time_function_gauss() should be modified as
described above (though it looks like that to me), because I am not much
familiar with adjoint or GPU simulations. So please look for all
comments containing "FS FS". Thanks.


There is one little thing though when using a point source and not a
Ricker wavelet:
The value of half duration is read from the f0 value in FORCESOLUTION,
actually meaning a value of half duration and not frequency. In my
opinion this should not be a problem, since in old SPECFEM3d versions
this confusion was similar when reading the value of center frequency
from the hdur line in file CMTSOLUTION.
Also, when setting f0 = 0.0 in FORCESOLUTION, the variable hdur is set
to TINYVAL and not to 5*DT which is done in the analogous case of using
a Moment tensor source with a Heaviside. But the user could simply write
f0=5*DT to FORCESOLUTION when he needs a steep gaussian.


I hope this could help. I am happy about any feedback in case there is a
better (more suited) solution to this. Thanks!

Florian

-- 
--------------------------------------------------------
 Dr. Florian Schumacher      <florian.schumacher at rub.de>
 Ruhr-Universität Bochum
 Institute of Geology, Mineralogy and Geophysics
 NA 3/168
 44780 Bochum
 Germany                      Tel. +49 (0)234 32 23273
--------------------------------------------------------
-------------- next part --------------
A non-text attachment was scrubbed...
Name: setup_sources_receivers.f90
Type: text/x-fortran
Size: 41179 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/cig-seismo/attachments/20150903/eee05265/attachment-0002.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: compute_add_sources_viscoelastic.f90
Type: text/x-fortran
Size: 42867 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/cig-seismo/attachments/20150903/eee05265/attachment-0003.bin>