[cig-commits] r20736 - in seismo/3D/SPECFEM3D/trunk: . doc/USER_MANUAL examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files examples/Mount_StHelens/in_data_files examples/homogeneous_halfspace/in_data_files examples/homogeneous_poroelastic/in_data_files examples/layered_halfspace examples/layered_halfspace/in_data_files examples/meshfem3D_examples/many_interfaces examples/meshfem3D_examples/simple_model/in_data_files examples/meshfem3D_examples/socal1D/in_data_files examples/tomographic_model/in_data_files examples/waterlayered_halfspace/in_data_files in_data_files src/generate_databases src/shared src/specfem3D
joseph.charles at geodynamics.org
joseph.charles at geodynamics.org
Wed Sep 19 11:47:13 PDT 2012
Author: joseph.charles
Date: 2012-09-19 11:47:12 -0700 (Wed, 19 Sep 2012)
New Revision: 20736
Modified:
seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.pdf
seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.tex
seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/homogeneous_poroelastic/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/README
seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/process.sh
seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/many_interfaces/Par_file
seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/socal1D/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90
seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90
seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90
seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90
seismo/3D/SPECFEM3D/trunk/src/shared/sum_kernels.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt
Log:
adds point source flags in Par_file related to the use of an inclined force
Modified: seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.tex
===================================================================
--- seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.tex 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.tex 2012-09-19 18:47:12 UTC (rev 20736)
@@ -1054,6 +1054,18 @@
than requested) seismograms in this directory. On a fast machine set
\texttt{NTSTEP\_BETWEEN\_OUTPUT\_SEISMOS} to a relatively high value
to avoid writing to the seismograms too often. This feature is only relevant for the solver.
+\item [{\texttt{USE\_FORCE\_POINT\_SOURCE}}] Turn this flag on to use a force source
+instead of a \texttt{CMTSOLUTION} moment-tensor source. This can be useful e.g. for oil industry
+foothills simulations in which the source is a vertical force, normal force,
+tilted force, or an impact etc. You will need to give the East, North and vertical components of the force vector;
+thus refer to Appendix~\ref{cha:Coordinates} for the orientation of the reference frame.
+\item [{\texttt{FACTOR\_FORCE\_SOURCE}}] This parameter specifies the magnitude of the force source.
+\item [{\texttt{COMPONENT\_DIR\_VECT\_SOURCE\_E}}] This parameter specifies the East component of
+an arbitrary (non-unitary) direction vector of the force source.
+\item [{\texttt{COMPONENT\_DIR\_VECT\_SOURCE\_N}}] This parameter specifies the North component of
+an arbitrary (non-unitary) direction vector of the force source.
+\item [{\texttt{COMPONENT\_DIR\_VECT\_SOURCE\_ZUP}}] This parameter specifies the vertical component of
+an arbitrary (non-unitary) direction vector of the force source.
\item [{\texttt{PRINT\_SOURCE\_TIME\_FUNCTION}}] Turn this flag on to print
information about the source time function in the file \texttt{in\_out\_files/OUTPUT\_FILES/plot\_source\_time\_function.txt}.
This feature is only relevant for the solver.
@@ -1188,7 +1200,8 @@
%
\begin{figure}
\noindent \begin{centering}
-\includegraphics[width=3in]{figures/gauss_vs_triangle_mod.eps}
+%%%\includegraphics[width=3in]{figures/gauss_vs_triangle_mod.eps}
+\includegraphics[width=3in]{figures/gauss_vs_triangle_mod.pdf}
\par\end{centering}
\caption{Comparison of the shape of a triangle and the Gaussian function actually
@@ -1231,7 +1244,7 @@
{\small }%
\begin{figure}[H]
\noindent \begin{centering}
-{\small \includegraphics[width=5in]{figures/source_timing.eps} }
+{\small \includegraphics[width=5in]{figures/source_timing.pdf} }
\par\end{centering}{\small \par}
\caption{Example of timing for three sources. The center of the first source
@@ -2848,7 +2861,7 @@
\begin{figure}[ht]
\noindent \begin{centering}
-\includegraphics[scale=0.6]{figures/IRIS_band_codes.eps}
+\includegraphics[scale=0.6]{figures/IRIS_band_codes.pdf}
\par\end{centering}
\caption{The band code convention is based on the sampling rate and the response band of instruments. Please visit \texttt{www.iris.edu/manuals/SEED\_appA.htm} for further information. Grey rows show the relative band-code range in SPECFEM3D, and the band codes used to name SEM seismograms are denoted in red.}
Modified: seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files/Par_file 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files/Par_file 2012-09-19 18:47:12 UTC (rev 20736)
@@ -63,11 +63,13 @@
# use a force source located exactly at a grid point instead of a CMTSOLUTION source
# this can be useful e.g. for oil industry foothills simulations or asteroid simulations
-# in which the source is a vertical force, normal force, impact etc.
+# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
FACTOR_FORCE_SOURCE = 1.d15
-# direction in comp E/N/Z = 1/2/3
-COMPONENT_FORCE_SOURCE = 3
+# components of a (non-unitary) direction vector for the force source on the E/N/Z_UP basis
+COMPONENT_DIR_VECT_SOURCE_E = 1.d0
+COMPONENT_DIR_VECT_SOURCE_N = -2.d0
+COMPONENT_DIR_VECT_SOURCE_Z_UP = -1.d0
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .true.
Modified: seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files/Par_file 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files/Par_file 2012-09-19 18:47:12 UTC (rev 20736)
@@ -63,11 +63,13 @@
# use a force source located exactly at a grid point instead of a CMTSOLUTION source
# this can be useful e.g. for oil industry foothills simulations or asteroid simulations
-# in which the source is a vertical force, normal force, impact etc.
+# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
FACTOR_FORCE_SOURCE = 1.d15
-# direction in comp E/N/Z = 1/2/3
-COMPONENT_FORCE_SOURCE = 3
+# components of a (non-unitary) direction vector for the force source on the E/N/Z_UP basis
+COMPONENT_DIR_VECT_SOURCE_E = 1.d0
+COMPONENT_DIR_VECT_SOURCE_N = -2.d0
+COMPONENT_DIR_VECT_SOURCE_Z_UP = -1.d0
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .true.
Modified: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/Par_file 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/Par_file 2012-09-19 18:47:12 UTC (rev 20736)
@@ -63,11 +63,13 @@
# use a force source located exactly at a grid point instead of a CMTSOLUTION source
# this can be useful e.g. for oil industry foothills simulations or asteroid simulations
-# in which the source is a vertical force, normal force, impact etc.
+# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
FACTOR_FORCE_SOURCE = 1.d15
-# direction in comp E/N/Z = 1/2/3
-COMPONENT_FORCE_SOURCE = 3
+# components of a (non-unitary) direction vector for the force source on the E/N/Z_UP basis
+COMPONENT_DIR_VECT_SOURCE_E = 1.d0
+COMPONENT_DIR_VECT_SOURCE_N = -2.d0
+COMPONENT_DIR_VECT_SOURCE_Z_UP = -1.d0
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .false.
Modified: seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/in_data_files/Par_file 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/in_data_files/Par_file 2012-09-19 18:47:12 UTC (rev 20736)
@@ -65,11 +65,13 @@
# use a force source located exactly at a grid point instead of a CMTSOLUTION source
# this can be useful e.g. for oil industry foothills simulations or asteroid simulations
-# in which the source is a vertical force, normal force, impact etc.
+# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
FACTOR_FORCE_SOURCE = 1.d15
-# direction in comp E/N/Z = 1/2/3
-COMPONENT_FORCE_SOURCE = 3
+# components of a (non-unitary) direction vector for the force source on the E/N/Z_UP basis
+COMPONENT_DIR_VECT_SOURCE_E = 1.d0
+COMPONENT_DIR_VECT_SOURCE_N = -2.d0
+COMPONENT_DIR_VECT_SOURCE_Z_UP = -1.d0
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .false.
Modified: seismo/3D/SPECFEM3D/trunk/examples/homogeneous_poroelastic/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/homogeneous_poroelastic/in_data_files/Par_file 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/homogeneous_poroelastic/in_data_files/Par_file 2012-09-19 18:47:12 UTC (rev 20736)
@@ -65,11 +65,13 @@
# use a force source located exactly at a grid point instead of a CMTSOLUTION source
# this can be useful e.g. for oil industry foothills simulations or asteroid simulations
-# in which the source is a vertical force, normal force, impact etc.
+# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
FACTOR_FORCE_SOURCE = 1.d15
-# direction in comp E/N/Z = 1/2/3
-COMPONENT_FORCE_SOURCE = 3
+# components of a (non-unitary) direction vector for the force source on the E/N/Z_UP basis
+COMPONENT_DIR_VECT_SOURCE_E = 1.d0
+COMPONENT_DIR_VECT_SOURCE_N = -2.d0
+COMPONENT_DIR_VECT_SOURCE_Z_UP = -1.d0
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .true.
Modified: seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/README
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/README 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/README 2012-09-19 18:47:12 UTC (rev 20736)
@@ -107,14 +107,6 @@
5. run simulation:
- - modify src/shared/constants.h:
- to use a vertical force source, with a Ricker wavelet source time function,
- turn flag on for parameter USE_FORCE_POINT_SOURCE:
- ...
- logical, parameter :: USE_FORCE_POINT_SOURCE = .true.
- double precision, parameter :: FACTOR_FORCE_SOURCE = -1.d15
- ...
-
- compile specfem3D:
> make xspecfem3D
Modified: seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/in_data_files/Par_file 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/in_data_files/Par_file 2012-09-19 18:47:12 UTC (rev 20736)
@@ -64,11 +64,13 @@
# use a force source located exactly at a grid point instead of a CMTSOLUTION source
# this can be useful e.g. for oil industry foothills simulations or asteroid simulations
-# in which the source is a vertical force, normal force, impact etc.
+# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
FACTOR_FORCE_SOURCE = 1.d15
-# direction in comp E/N/Z = 1/2/3
-COMPONENT_FORCE_SOURCE = 3
+# components of a (non-unitary) direction vector for the force source on the E/N/Z_UP basis
+COMPONENT_DIR_VECT_SOURCE_E = 1.d0
+COMPONENT_DIR_VECT_SOURCE_N = -2.d0
+COMPONENT_DIR_VECT_SOURCE_Z_UP = -1.d0
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .false.
Modified: seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/process.sh 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/process.sh 2012-09-19 18:47:12 UTC (rev 20736)
@@ -34,22 +34,6 @@
rm -rf in_out_files/OUTPUT_FILES/*
rm -rf in_out_files/DATABASES_MPI/*
-# compiles executables in root directory
-cd ../../
-
-# compiles with flag for a point force (with a ricker source time function)
-cd src/shared/
-cp constants.h constants.h.org
-sed -e "s:USE_FORCE_POINT_SOURCE.*:USE_FORCE_POINT_SOURCE = .true. :" constants.h.org > constants.h.tmp
-sed -e "s:FACTOR_FORCE_SOURCE.*:FACTOR_FORCE_SOURCE = -1.d15 :" constants.h.tmp > constants.h
-rm -f constants.h.tmp
-cd ../../
-make > tmp.log
-
-# backup & restores original file again
-cp src/shared/constants.h $currentdir/in_out_files/OUTPUT_FILES/constants.h.bak
-mv src/shared/constants.h.org src/shared/constants.h
-
cd $currentdir
# links executables
Modified: seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/many_interfaces/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/many_interfaces/Par_file 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/many_interfaces/Par_file 2012-09-19 18:47:12 UTC (rev 20736)
@@ -64,11 +64,13 @@
# use a force source located exactly at a grid point instead of a CMTSOLUTION source
# this can be useful e.g. for oil industry foothills simulations or asteroid simulations
-# in which the source is a vertical force, normal force, impact etc.
+# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
FACTOR_FORCE_SOURCE = 1.d15
-# direction in comp E/N/Z = 1/2/3
-COMPONENT_FORCE_SOURCE = 3
+# components of a (non-unitary) direction vector for the force source on the E/N/Z_UP basis
+COMPONENT_DIR_VECT_SOURCE_E = 1.d0
+COMPONENT_DIR_VECT_SOURCE_N = -2.d0
+COMPONENT_DIR_VECT_SOURCE_Z_UP = -1.d0
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .false.
Modified: seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/in_data_files/Par_file 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/in_data_files/Par_file 2012-09-19 18:47:12 UTC (rev 20736)
@@ -64,11 +64,13 @@
# use a force source located exactly at a grid point instead of a CMTSOLUTION source
# this can be useful e.g. for oil industry foothills simulations or asteroid simulations
-# in which the source is a vertical force, normal force, impact etc.
+# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
FACTOR_FORCE_SOURCE = 1.d15
-# direction in comp E/N/Z = 1/2/3
-COMPONENT_FORCE_SOURCE = 3
+# components of a (non-unitary) direction vector for the force source on the E/N/Z_UP basis
+COMPONENT_DIR_VECT_SOURCE_E = 1.d0
+COMPONENT_DIR_VECT_SOURCE_N = -2.d0
+COMPONENT_DIR_VECT_SOURCE_Z_UP = -1.d0
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .false.
Modified: seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/socal1D/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/socal1D/in_data_files/Par_file 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/socal1D/in_data_files/Par_file 2012-09-19 18:47:12 UTC (rev 20736)
@@ -64,11 +64,13 @@
# use a force source located exactly at a grid point instead of a CMTSOLUTION source
# this can be useful e.g. for oil industry foothills simulations or asteroid simulations
-# in which the source is a vertical force, normal force, impact etc.
+# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
FACTOR_FORCE_SOURCE = 1.d15
-# direction in comp E/N/Z = 1/2/3
-COMPONENT_FORCE_SOURCE = 3
+# components of a (non-unitary) direction vector for the force source on the E/N/Z_UP basis
+COMPONENT_DIR_VECT_SOURCE_E = 1.d0
+COMPONENT_DIR_VECT_SOURCE_N = -2.d0
+COMPONENT_DIR_VECT_SOURCE_Z_UP = -1.d0
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .true.
Modified: seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/in_data_files/Par_file 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/in_data_files/Par_file 2012-09-19 18:47:12 UTC (rev 20736)
@@ -64,11 +64,13 @@
# use a force source located exactly at a grid point instead of a CMTSOLUTION source
# this can be useful e.g. for oil industry foothills simulations or asteroid simulations
-# in which the source is a vertical force, normal force, impact etc.
+# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
FACTOR_FORCE_SOURCE = 1.d15
-# direction in comp E/N/Z = 1/2/3
-COMPONENT_FORCE_SOURCE = 3
+# components of a (non-unitary) direction vector for the force source on the E/N/Z_UP basis
+COMPONENT_DIR_VECT_SOURCE_E = 1.d0
+COMPONENT_DIR_VECT_SOURCE_N = -2.d0
+COMPONENT_DIR_VECT_SOURCE_Z_UP = -1.d0
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .false.
Modified: seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/in_data_files/Par_file 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/in_data_files/Par_file 2012-09-19 18:47:12 UTC (rev 20736)
@@ -64,11 +64,13 @@
# use a force source located exactly at a grid point instead of a CMTSOLUTION source
# this can be useful e.g. for oil industry foothills simulations or asteroid simulations
-# in which the source is a vertical force, normal force, impact etc.
+# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
FACTOR_FORCE_SOURCE = 1.d15
-# direction in comp E/N/Z = 1/2/3
-COMPONENT_FORCE_SOURCE = 3
+# components of a (non-unitary) direction vector for the force source on the E/N/Z_UP basis
+COMPONENT_DIR_VECT_SOURCE_E = 1.d0
+COMPONENT_DIR_VECT_SOURCE_N = -2.d0
+COMPONENT_DIR_VECT_SOURCE_Z_UP = -1.d0
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .false.
Modified: seismo/3D/SPECFEM3D/trunk/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/in_data_files/Par_file 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/in_data_files/Par_file 2012-09-19 18:47:12 UTC (rev 20736)
@@ -65,11 +65,13 @@
# use a force source located exactly at a grid point instead of a CMTSOLUTION source
# this can be useful e.g. for oil industry foothills simulations or asteroid simulations
-# in which the source is a vertical force, normal force, impact etc.
+# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
FACTOR_FORCE_SOURCE = 1.d15
-# direction in comp E/N/Z = 1/2/3
-COMPONENT_FORCE_SOURCE = 3
+# components of a (non-unitary) direction vector for the force source on the E/N/Z_UP basis
+COMPONENT_DIR_VECT_SOURCE_E = 1.d0
+COMPONENT_DIR_VECT_SOURCE_N = -2.d0
+COMPONENT_DIR_VECT_SOURCE_Z_UP = -1.d0
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .false.
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -263,7 +263,8 @@
NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,FACTOR_FORCE_SOURCE, &
- COMPONENT_FORCE_SOURCE,IMODEL)
+ COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N, &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP,IMODEL)
! check that the code is running with the requested nb of processes
if(sizeprocs /= NPROC) then
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -56,9 +56,11 @@
! parameters read from parameter file
integer :: NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,SIMULATION_TYPE
- integer :: NSOURCES,COMPONENT_FORCE_SOURCE
+ integer :: NSOURCES
double precision :: DT,HDUR_MOVIE,FACTOR_FORCE_SOURCE
+ double precision :: COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N, &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP
logical :: ATTENUATION,USE_OLSEN_ATTENUATION, &
OCEANS,TOPOGRAPHY,SAVE_FORWARD,USE_FORCE_POINT_SOURCE
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -90,6 +90,8 @@
! for read_parameter_files
double precision :: DT
double precision :: HDUR_MOVIE,FACTOR_FORCE_SOURCE
+ double precision :: COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N, &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP
integer :: NPROC,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP, &
UTM_PROJECTION_ZONE,SIMULATION_TYPE
integer :: NSOURCES,NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY
@@ -101,7 +103,7 @@
logical :: ABSORBING_CONDITIONS,SAVE_FORWARD
logical :: ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
character(len=256) LOCAL_PATH
- integer :: COMPONENT_FORCE_SOURCE,IMODEL
+ integer :: IMODEL
! checks given arguments
print *
@@ -183,7 +185,8 @@
NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,FACTOR_FORCE_SOURCE, &
- COMPONENT_FORCE_SOURCE,IMODEL)
+ COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N, &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP,IMODEL)
print *, 'Slice list: '
print *, node_list(1:num_node)
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -91,6 +91,8 @@
integer NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
double precision DT
double precision HDUR_MOVIE,FACTOR_FORCE_SOURCE
+ double precision COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N, &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP
logical ATTENUATION,USE_OLSEN_ATTENUATION, &
OCEANS,TOPOGRAPHY,USE_FORCE_POINT_SOURCE
logical ABSORBING_CONDITIONS,SAVE_FORWARD
@@ -98,7 +100,7 @@
character(len=256) OUTPUT_FILES,LOCAL_PATH
integer NPROC
integer ier
- integer COMPONENT_FORCE_SOURCE,IMODEL
+ integer IMODEL
!--------------------------------------------
!!!! NL NL for external meshes
@@ -137,7 +139,8 @@
NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,FACTOR_FORCE_SOURCE, &
- COMPONENT_FORCE_SOURCE,IMODEL)
+ COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N, &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP,IMODEL)
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', OUTPUT_FILES_PATH(1:len_trim(OUTPUT_FILES_PATH)))
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -32,7 +32,8 @@
NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION,NTSTEP_BETWEEN_OUTPUT_INFO, &
SIMULATION_TYPE,SAVE_FORWARD,NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
- USE_FORCE_POINT_SOURCE,FACTOR_FORCE_SOURCE,COMPONENT_FORCE_SOURCE,IMODEL )
+ USE_FORCE_POINT_SOURCE,FACTOR_FORCE_SOURCE,COMPONENT_DIR_VECT_SOURCE_E, &
+ COMPONENT_DIR_VECT_SOURCE_N,COMPONENT_DIR_VECT_SOURCE_Z_UP,IMODEL)
implicit none
@@ -40,10 +41,12 @@
integer NPROC,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,SIMULATION_TYPE, NTSTEP_BETWEEN_READ_ADJSRC
integer NSOURCES,NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO,UTM_PROJECTION_ZONE
- integer NOISE_TOMOGRAPHY,COMPONENT_FORCE_SOURCE
+ integer NOISE_TOMOGRAPHY
integer IMODEL
double precision DT,HDUR_MOVIE,FACTOR_FORCE_SOURCE
+ double precision COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N
+ double precision COMPONENT_DIR_VECT_SOURCE_Z_UP
logical ATTENUATION,USE_OLSEN_ATTENUATION,OCEANS,TOPOGRAPHY,ABSORBING_CONDITIONS,SAVE_FORWARD
logical MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT,USE_HIGHRES_FOR_MOVIES
@@ -143,8 +146,12 @@
if(err_occurred() /= 0) return
call read_value_double_precision(FACTOR_FORCE_SOURCE, 'solver.FACTOR_FORCE_SOURCE')
if(err_occurred() /= 0) return
- call read_value_integer(COMPONENT_FORCE_SOURCE, 'solver.COMPONENT_FORCE_SOURCE')
+ call read_value_double_precision(COMPONENT_DIR_VECT_SOURCE_E, 'solver.COMPONENT_DIR_VECT_SOURCE_E')
if(err_occurred() /= 0) return
+ call read_value_double_precision(COMPONENT_DIR_VECT_SOURCE_N, 'solver.COMPONENT_DIR_VECT_SOURCE_N')
+ if(err_occurred() /= 0) return
+ call read_value_double_precision(COMPONENT_DIR_VECT_SOURCE_Z_UP, 'solver.COMPONENT_DIR_VECT_SOURCE_Z_UP')
+ if(err_occurred() /= 0) return
call read_value_logical(PRINT_SOURCE_TIME_FUNCTION, 'solver.PRINT_SOURCE_TIME_FUNCTION')
if(err_occurred() /= 0) return
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -98,6 +98,8 @@
! for read_parameter_files
double precision :: DT
double precision :: HDUR_MOVIE,FACTOR_FORCE_SOURCE
+ double precision :: COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N, &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP
integer :: NPROC,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP, &
UTM_PROJECTION_ZONE,SIMULATION_TYPE
integer :: NSOURCES,NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY
@@ -109,7 +111,7 @@
logical :: ABSORBING_CONDITIONS,SAVE_FORWARD
logical :: ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
character(len=256) LOCAL_PATH
- integer :: COMPONENT_FORCE_SOURCE,IMODEL
+ integer :: IMODEL
! smoothing parameters
character(len=256) :: ks_file
@@ -221,7 +223,8 @@
NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,FACTOR_FORCE_SOURCE, &
- COMPONENT_FORCE_SOURCE,IMODEL)
+ COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N, &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP,IMODEL)
! checks if number of MPI process as specified
if (sizeprocs /= NPROC) then
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/sum_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/sum_kernels.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/sum_kernels.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -104,11 +104,13 @@
! for read_parameter_files
double precision :: DT
double precision :: HDUR_MOVIE,FACTOR_FORCE_SOURCE
+ double precision :: COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N, &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP
integer :: NPROC,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP, &
UTM_PROJECTION_ZONE,SIMULATION_TYPE
integer :: NSOURCES,NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY
integer :: NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
- integer :: COMPONENT_FORCE_SOURCE,IMODEL
+ integer :: IMODEL
logical :: MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
USE_HIGHRES_FOR_MOVIES,SUPPRESS_UTM_PROJECTION
logical :: ATTENUATION,USE_OLSEN_ATTENUATION, &
@@ -160,7 +162,8 @@
NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,FACTOR_FORCE_SOURCE, &
- COMPONENT_FORCE_SOURCE,IMODEL)
+ COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N, &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP,IMODEL)
! checks if number of MPI process as specified
if (sizeprocs /= NPROC) then
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -140,7 +140,7 @@
! write(*,*) "fortran dt = ", dt
! change dt -> DT
call compute_add_sources_ac_cuda(Mesh_pointer, phase_is_inner, &
- NSOURCES,stf_pre_compute, myrank)
+ NSOURCES, stf_pre_compute, myrank)
endif
else ! .NOT. GPU_MODE
@@ -432,7 +432,6 @@
endif
enddo
stf_used_total = stf_used_total + sum(stf_pre_compute(:))
-
! only implements SIMTYPE=3
call compute_add_sources_ac_s3_cuda(Mesh_pointer, phase_is_inner, &
NSOURCES,stf_pre_compute, myrank)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -29,7 +29,7 @@
subroutine compute_add_sources_elastic( NSPEC_AB,NGLOB_AB,accel, &
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
- xi_source,eta_source,gamma_source,nu_source, &
+ xi_source,eta_source,gamma_source, &
hdur,hdur_gaussian,tshift_cmt,dt,t0,sourcearrays, &
ispec_is_elastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
@@ -46,7 +46,7 @@
mask_noise,noise_surface_movie, &
nrec_local,number_receiver_global, &
nsources_local,USE_FORCE_POINT_SOURCE, &
- FACTOR_FORCE_SOURCE,COMPONENT_FORCE_SOURCE
+ FACTOR_FORCE_SOURCE
implicit none
@@ -68,7 +68,6 @@
integer :: NSOURCES,myrank,it
integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
- double precision, dimension(3,3,NSOURCES) :: nu_source
double precision, dimension(NSOURCES) :: hdur,hdur_gaussian,tshift_cmt
double precision :: dt,t0
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrays
@@ -122,21 +121,29 @@
if (SIMULATION_TYPE == 1 .and. NOISE_TOMOGRAPHY == 0 .and. nsources_local > 0) then
if(GPU_MODE) then
- do isource = 1,NSOURCES
- if( USE_RICKER_IPATI ) then
- stf_pre_compute(isource) = comp_source_time_function_rickr( &
- dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
- else
- stf_pre_compute(isource) = comp_source_time_function( &
- dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
- endif
- enddo
- ! only implements SIMTYPE=1 and NOISE_TOM=0
- ! write(*,*) "fortran dt = ", dt
- ! change dt -> DT
- call compute_add_sources_el_cuda(Mesh_pointer, &
- phase_is_inner,NSOURCES, &
- stf_pre_compute, myrank)
+ if( NSOURCES > 0 ) then
+ do isource = 1,NSOURCES
+ if(USE_FORCE_POINT_SOURCE) then
+ ! precomputes source time function factor
+ stf_pre_compute(isource) = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr( &
+ dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
+ else
+ if( USE_RICKER_IPATI ) then
+ stf_pre_compute(isource) = comp_source_time_function_rickr( &
+ dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
+ else
+ stf_pre_compute(isource) = comp_source_time_function( &
+ dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ endif
+ endif
+ enddo
+ ! only implements SIMTYPE=1 and NOISE_TOM=0
+ ! write(*,*) "fortran dt = ", dt
+ ! change dt -> DT
+ call compute_add_sources_el_cuda(Mesh_pointer, phase_is_inner, &
+ NSOURCES, stf_pre_compute, myrank)
+
+ endif
else ! .NOT. GPU_MODE
@@ -168,15 +175,25 @@
!endif
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = FACTOR_FORCE_SOURCE * &
- comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+ stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
- ! we use a force in a single direction along one of the components:
- ! x/y/z or E/N/Z-direction would correspond to 1/2/3 = COMPONENT_FORCE_SOURCE
- ! e.g. nu_source(:,3) here would be a source normal to the surface (z-direction).
- accel(:,iglob) = accel(:,iglob) &
- + sngl( nu_source(COMPONENT_FORCE_SOURCE,:,isource) ) * stf_used
+ ! add the inclined force source array
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ accel(:,iglob) = accel(:,iglob) + sourcearrays(isource,:,i,j,k) * stf_used
+ enddo
+ enddo
+ enddo
+
else
if( USE_RICKER_IPATI) then
@@ -390,20 +407,29 @@
! adjoint simulations
if (SIMULATION_TYPE == 3 .and. NOISE_TOMOGRAPHY == 0 .and. nsources_local > 0) then
- if(GPU_MODE) then
- do isource = 1,NSOURCES
- if( USE_RICKER_IPATI ) then
- stf_pre_compute(isource) = comp_source_time_function_rickr( &
- dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
- else
- stf_pre_compute(isource) = comp_source_time_function( &
- dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ if(GPU_MODE) then
+ if( NSOURCES > 0 ) then
+ do isource = 1,NSOURCES
+ if(USE_FORCE_POINT_SOURCE) then
+ ! precomputes source time function factors
+ stf_pre_compute(isource) = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr( &
+ dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
+ else
+ if( USE_RICKER_IPATI ) then
+ stf_pre_compute(isource) = comp_source_time_function_rickr( &
+ dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
+ else
+ stf_pre_compute(isource) = comp_source_time_function( &
+ dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ endif
+ endif
+ enddo
+ ! only implements SIMTYPE=3
+ call compute_add_sources_el_s3_cuda(Mesh_pointer,stf_pre_compute, &
+ NSOURCES,phase_is_inner,myrank)
+
endif
- enddo
- call compute_add_sources_el_s3_cuda(Mesh_pointer,stf_pre_compute, &
- NSOURCES,phase_is_inner,myrank)
-
else ! .NOT. GPU_MODE
! backward source reconstruction
@@ -435,14 +461,25 @@
!endif
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = FACTOR_FORCE_SOURCE * &
- comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
+ stf = comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
- ! e.g. we use nu_source(:,3) here if we want a source normal to the surface.
- ! note: time step is now at NSTEP-it
- b_accel(:,iglob) = b_accel(:,iglob) &
- + sngl( nu_source(COMPONENT_FORCE_SOURCE,:,isource) ) * stf_used
+ ! add the inclined force source array
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ b_accel(:,iglob) = b_accel(:,iglob) + sourcearrays(isource,:,i,j,k) * stf_used
+ enddo
+ enddo
+ enddo
+
else
! see note above: time step corresponds now to NSTEP-it
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -32,7 +32,7 @@
rhoarraystore,phistore,tortstore,&
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
- xi_source,eta_source,gamma_source,nu_source, &
+ xi_source,eta_source,gamma_source, &
hdur,hdur_gaussian,tshift_cmt,dt,t0,sourcearrays, &
ispec_is_poroelastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
@@ -42,7 +42,8 @@
use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total, &
xigll,yigll,zigll,xi_receiver,eta_receiver,gamma_receiver,&
station_name,network_name,adj_source_file, &
- USE_FORCE_POINT_SOURCE,FACTOR_FORCE_SOURCE,COMPONENT_FORCE_SOURCE
+ USE_FORCE_POINT_SOURCE,FACTOR_FORCE_SOURCE
+
implicit none
include "constants.h"
@@ -66,7 +67,6 @@
integer :: NSOURCES,myrank,it
integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
- double precision, dimension(3,3,NSOURCES) :: nu_source
double precision, dimension(NSOURCES) :: hdur,hdur_gaussian,tshift_cmt
double precision :: dt,t0
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrays
@@ -154,19 +154,32 @@
endif
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+ stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
!stf_used = comp_source_time_function_gauss(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
- ! we use a force in a single direction along one of the components:
- ! x/y/z or E/N/Z-direction would correspond to 1/2/3 = COMPONENT_FORCE_SOURCE
- ! e.g. nu_source(:,3) here would be a source normal to the surface (z-direction).
+ ! add the inclined force source array
! the source is applied to both solid and fluid phase: bulk source.
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
! solid phase
- accels(:,iglob) = accels(:,iglob) &
- + (1._CUSTOM_REAL - phil/tortl) * sngl( nu_source(COMPONENT_FORCE_SOURCE,:,isource) ) * stf_used
+ accels(:,iglob) = accels(:,iglob) + &
+ (1._CUSTOM_REAL - phil/tortl) * sourcearrays(isource,:,i,j,k) * stf_used
! fluid phase
- accelw(:,iglob) = accelw(:,iglob) &
- + (1._CUSTOM_REAL - rhol_f/rhol_bar) * sngl( nu_source(COMPONENT_FORCE_SOURCE,:,isource) ) * stf_used
+ accelw(:,iglob) = accelw(:,iglob) + &
+ (1._CUSTOM_REAL - rhol_f/rhol_bar) * sourcearrays(isource,:,i,j,k) * stf_used
+ enddo
+ enddo
+ enddo
else
@@ -402,17 +415,32 @@
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
! should be the same than for the forward simulation (check above)
- stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
+ stf = comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
- ! e.g. we use nu_source(:,3) here if we want a source normal to the surface.
+ ! add the inclined force source array
+ ! the source is applied to both solid and fluid phase: bulk source
! note: time step is now at NSTEP-it
- ! the source is applied to both solid and fluid phase: bulk source
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec_selected_source(isource))
! solid phase
- b_accels(:,iglob) = b_accels(:,iglob) &
- + (1._CUSTOM_REAL - phil/tortl) * sngl( nu_source(COMPONENT_FORCE_SOURCE,:,isource) ) * stf_used
+ b_accels(:,iglob) = b_accels(:,iglob) + &
+ (1._CUSTOM_REAL - phil/tortl) * sourcearrays(isource,:,i,j,k) * stf_used
! fluid phase
- b_accelw(:,iglob) = b_accelw(:,iglob) &
- + (1._CUSTOM_REAL - rhol_f/rhol_bar) * sngl( nu_source(COMPONENT_FORCE_SOURCE,:,isource) ) * stf_used
+ b_accelw(:,iglob) = b_accelw(:,iglob) + &
+ (1._CUSTOM_REAL - rhol_f/rhol_bar) * sourcearrays(isource,:,i,j,k) * stf_used
+ enddo
+ enddo
+ enddo
else
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.F90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.F90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -235,7 +235,7 @@
call compute_add_sources_elastic( NSPEC_AB,NGLOB_AB,accel, &
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
- xi_source,eta_source,gamma_source,nu_source, &
+ xi_source,eta_source,gamma_source, &
hdur,hdur_gaussian,tshift_cmt,dt,t0,sourcearrays, &
ispec_is_elastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -205,7 +205,7 @@
rhoarraystore,phistore,tortstore,&
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
- xi_source,eta_source,gamma_source,nu_source, &
+ xi_source,eta_source,gamma_source, &
hdur,hdur_gaussian,tshift_cmt,dt,t0,sourcearrays, &
ispec_is_poroelastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -48,7 +48,8 @@
NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,FACTOR_FORCE_SOURCE, &
- COMPONENT_FORCE_SOURCE,IMODEL)
+ COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N, &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP,IMODEL)
! GPU_MODE is in par_file
call read_gpu_mode(GPU_MODE,GRAVITY)
@@ -265,6 +266,13 @@
stop 'ABSORBING_CONDITIONS must have NGLLX = NGLLY = NGLLZ'
endif
+ ! inclined force source
+ if( USE_FORCE_POINT_SOURCE ) then
+ if( COMPONENT_DIR_VECT_SOURCE_E .eq. 0.d0 .and. COMPONENT_DIR_VECT_SOURCE_N .eq. 0.d0 &
+ .and. COMPONENT_DIR_VECT_SOURCE_Z_UP .eq. 0.d0 ) &
+ stop 'USE_FORCE_POINT_SOURCE requires a non null direction vector'
+ endif
+
! exclusive movie flags
if( EXTERNAL_MESH_MOVIE_SURFACE .or. EXTERNAL_MESH_CREATE_SHAKEMAP ) then
if( EXTERNAL_MESH_MOVIE_SURFACE .and. EXTERNAL_MESH_CREATE_SHAKEMAP ) &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -39,18 +39,21 @@
nu_source,iglob_is_surface_external_mesh,ispec_is_surface_external_mesh, &
ispec_is_acoustic,ispec_is_elastic,ispec_is_poroelastic, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk, &
- USE_FORCE_POINT_SOURCE,COMPONENT_FORCE_SOURCE)
+ USE_FORCE_POINT_SOURCE,COMPONENT_DIR_VECT_SOURCE_E, &
+ COMPONENT_DIR_VECT_SOURCE_N,COMPONENT_DIR_VECT_SOURCE_Z_UP)
implicit none
include "constants.h"
integer NPROC,UTM_PROJECTION_ZONE
- integer NSPEC_AB,NGLOB_AB,NSOURCES,COMPONENT_FORCE_SOURCE
+ integer NSPEC_AB,NGLOB_AB,NSOURCES
logical PRINT_SOURCE_TIME_FUNCTION,SUPPRESS_UTM_PROJECTION,USE_FORCE_POINT_SOURCE
double precision DT
+ double precision COMPONENT_DIR_VECT_SOURCE_E,COMPONENT_DIR_VECT_SOURCE_N, &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
@@ -729,7 +732,9 @@
write(IMAIN,*) ' j index of source in that element: ',nint(eta_source(isource))
write(IMAIN,*) ' k index of source in that element: ',nint(gamma_source(isource))
write(IMAIN,*)
- write(IMAIN,*) ' component direction: ',COMPONENT_FORCE_SOURCE
+ write(IMAIN,*) ' component of direction vector in East direction: ',COMPONENT_DIR_VECT_SOURCE_E
+ write(IMAIN,*) ' component of direction vector in North direction: ',COMPONENT_DIR_VECT_SOURCE_N
+ write(IMAIN,*) ' component of direction vector in Vertical direction: ',COMPONENT_DIR_VECT_SOURCE_Z_UP
write(IMAIN,*)
write(IMAIN,*) ' nu1 = ',nu_source(1,:,isource)
write(IMAIN,*) ' nu2 = ',nu_source(2,:,isource)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -117,7 +117,8 @@
nu_source,iglob_is_surface_external_mesh,ispec_is_surface_external_mesh,&
ispec_is_acoustic,ispec_is_elastic,ispec_is_poroelastic, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk, &
- USE_FORCE_POINT_SOURCE,COMPONENT_FORCE_SOURCE)
+ USE_FORCE_POINT_SOURCE,COMPONENT_DIR_VECT_SOURCE_E, &
+ COMPONENT_DIR_VECT_SOURCE_N,COMPONENT_DIR_VECT_SOURCE_Z_UP)
if(abs(minval(tshift_cmt)) > TINYVAL) call exit_MPI(myrank,'one tshift_cmt must be zero, others must be positive')
@@ -600,7 +601,14 @@
endif
! elastic source
if( ispec_is_elastic(ispec) ) then
- sourcearray(:,i,j,k) = nu_source(COMPONENT_FORCE_SOURCE,:,isource)
+ ! we use an inclined force defined by its magnitude and the projections
+ ! of an arbitrary (non-unitary) direction vector on the E/N/Z_UP basis:
+ sourcearray(:,i,j,k) = FACTOR_FORCE_SOURCE * &
+ ( nu_source(1,:,isource) * COMPONENT_DIR_VECT_SOURCE_E + &
+ nu_source(2,:,isource) * COMPONENT_DIR_VECT_SOURCE_N + &
+ nu_source(3,:,isource) * COMPONENT_DIR_VECT_SOURCE_Z_UP ) / &
+ sqrt( COMPONENT_DIR_VECT_SOURCE_E**2 + COMPONENT_DIR_VECT_SOURCE_N**2 + &
+ COMPONENT_DIR_VECT_SOURCE_Z_UP**2 )
endif
endif
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-09-19 18:47:12 UTC (rev 20736)
@@ -157,7 +157,9 @@
! parameters for a force source located exactly at a grid point
logical :: USE_FORCE_POINT_SOURCE
double precision :: FACTOR_FORCE_SOURCE
- integer :: COMPONENT_FORCE_SOURCE
+ double precision :: COMPONENT_DIR_VECT_SOURCE_E
+ double precision :: COMPONENT_DIR_VECT_SOURCE_N
+ double precision :: COMPONENT_DIR_VECT_SOURCE_Z_UP
! parameters
integer :: SIMULATION_TYPE
Modified: seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt 2012-09-19 17:23:39 UTC (rev 20735)
+++ seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt 2012-09-19 18:47:12 UTC (rev 20736)
@@ -98,240 +98,6 @@
easy to do once we have an expression for the weighting factors
------------------------------------------------
-Par_file modifications:
-------------------------------------------------
-
-- suggestion 05:
-----------------
-
-Subject: Re: [CIG-SEISMO] simulation of seismic wave propagation excited by a force
-Date: Tue, 11 Sep 2012 18:34:21 +0200
-From: Dimitri Komatitsch
-To: Christina Morency
-CC: Yang Luo, Yingzi Ying
-
-Hi Christina and Yang, Hi all,
-
-OK, great. This is excellent news.
-
-So let us just move this "USE_FORCE_POINT_SOURCE" flag from
-setup/constants.h.in to DATA/Par_file, and let us mention it in the manual.
-
-Jo, could you please do that and commit the changes?
-
-Thank you,
-Dimitri.
-
-On 09/11/2012 06:17 PM, Christina Morency wrote:
-> Hi Dimitri,
->
-> Yes, I can confirm what Yang is saying, it is already in the package,
-> with interpolation if the source is not located on a GLL point (as it is
-> indeed done in the 2D btw).
->
-> Christina
->
-> On 09/11/2012 09:14 AM, Yang Luo wrote:
->> I think Daniel (or someone else) has already copied necessary lines from old package.
->>
->> There is an option "USE_FORCE_POINT_SOURCE" in constant.h, which allows one to use a point force
->> (which Yingzi already figured out, or maybe it's already in the manual)
->> logical, parameter :: USE_FORCE_POINT_SOURCE = .false.
->> double precision, parameter :: FACTOR_FORCE_SOURCE = 1.d15
->> integer, parameter :: COMPONENT_FORCE_SOURCE = 3 ! takes direction in comp E/N/Z = 1/2/3
->>
->> Interpolation is included, if the point force is not located exactly on a GLL point.
-
-=================
-
-add an option to use a point force source in SPECFEM3D instead of a CMTSOLUTION; some users need that (very classical for instance
-in non destructive testing) and it is easy to do; the only thing to decide is what convention to use in DATA/Par_file: add a flag
-for that?
-put a flag to choose between a Ricker source and a Heaviside source when a force source is used
-(we can cut and paste the Heaviside implementation from the case of a CMT source)
-
-Feedback from Qinya: I think this is a great idea. In the past, we always had multiple versions of the same specfem code, just to be able to run point force sources. We can add a parameter `SOURCE_TYPE' in Par file, which is = 0 for moment-tensor source ( reads CMTSOLUTION, and assumes Heavside/err as source time function), and = 1 for point force source (and read something like PNTSOLUTION, including force location, direction, time shift and half duration). Now there is no standard stf to use for point force, but Heavside function does not seem to be a good idea because it will produce a static offset on all seismograms, which is not desirable. On the other hand, I don't see why Richer wavelet (2nd order derivative of Gaussian) is necessary better than Gaussian itself. One way to give user the freedom of choice is to put the stf info also in the PNTSOLUTION.
-
-Subject: Re: [CIG-SEISMO] simulation of seismic wave propagation excited by a force
-Date: Tue, 11 Sep 2012 18:01:38 +0200
-From: Dimitri Komatitsch
-To: Yingzi Ying
-
-Dear Yingzi,
-
-You are right, in the current code there is no option in the Par_file to
-use a point source, but we added this to the todo list right before the
-summer and thus we should/will definitely do it. It is easy to do: just
-add the Ricker wavelet to the accel() vector (which is in fact a force
-vector before we divide it by the mass matrix), to the component you
-want (vertical or horizontal) and you are all set.
-
-Joseph, could you maybe do it? below are a few lines from one of my test
-codes, in which I do just that; let us cut and paste these lines in
-SPECFEM3D (it is just a matter of selecting the right grid point; which
-we already do for a CMT source, thus we could add a force instead of a
-CMT at that location; the only difference being that the CMT source is
-added to all the points of a given spectral element, while a force
-vector should be added to a single point).
-
-Let me talk to Joseph tomorrow to see if we can release a patch of the
-source code soon.
-
-I also cc Yang to see if he has already done that in a local version at
-Princeton maybe (for instance to model a force for some of his oil
-industry simulations in foothill regions?).
-
-Another difficulty is if the location of the force does not correspond
-to a grid point; then there are two options:
-
-- select the closest grid point; i.e. slightly change the location of
-the source (can be done automatically by the code)
-
-- use Christina Morency's way of adding a force source between grid
-points; I think I remember she did that in the 2D code, but I am not
-sure how; thus I cc her.
-
-(that second option would be better, if Christina can confirm that it
-works fine)
-
-Best wishes,
-Dimitri.
-
-==================================================================
-
-! ipoinsource is the global grid point at which we add the source
-! "3" means a vertical force (i.e. we add the force to the vertical
-! component of the force vector, which is later divided by the mass
-! matrix to get the acceleration vector)
-! factor_force is the norm of the force vector to add
-! ricker(t) is the Ricker wavelet at time t = (it-1)*deltat
-
- accel(3,ipoinsource) = accel(3,ipoinsource) - factor_force * ricker(t)
-
-==================================================================
-
-On 09/11/2012 12:59 PM, Yingzi Ying wrote:
-> Hi all,
->
-> I am trying to use specfem3d to simulate waves excited by a vertical
-> force with user-defined source forms such as ricker wavelet.
->
-> I tried to do it through adjoint simulation by setting 'SIMULATION_TYPE
-> = 2' and put source position at STATIONS_ADJOINT and receivers at
-> CMTSOLUTION. But I failed.
->
-> Does anyone have an idea or suggestion how to do active simulations with
-> specfem3d?
->
-> Thanks in advance.
->
-> Best regards,
-> Yingzi
->
-
-
--------- Original Message --------
-Subject: Re: force in SPECFEM3D
-Date: Wed, 18 Jul 2012 23:50:19 +0200
-From: Dimitri Komatitsch
-To: Carl Tape, daniel peter
-
-Hi all,
-
-Yes, let us definitely add it and make it an option in the Par_file.
-This will be very useful also in the case of non destructive testing of
-materials (in which case the source is almost always a set of forces).
-
-What probably exists is a routine to add a force at a GLL grid point
-(which is easy to do). If we also want to have an option to model a
-source located somewhere between grid points, then I guess we will need
-to see how the discrete source term looks like (I think Daniel once told
-me that Christina had done it in the 2D code; probably summing the
-Lagrange interpolants at all the GLL points of the element? (as we do
-for receivers located between grid points?).
-
-Thanks,
-Dimitri.
-
-On 07/18/2012 11:01 PM, Jeroen Tromp wrote:
-> Hi Daniel:
->
-> The adjoint source involves a "force", so the necessary code exists. But
-> it can't hurt to explicitly allow a source that is a force rather than a
-> moment tensor.
->
-> Best regards,
->
-> Jeroen
->
-> On 7/18/12 8:56 PM, Carl Tape wrote:
->> Hi Daniel, Dimitri, Jeroen--
->>
->> I am continuing some benchmark simulations for static offsets in
->> SPECFEM3D: 0-D, 1-D, 3-D models. One of the applications is to make a
->> Green's function database similar to what is used in GPS inversions of
->> fault slip. This means using a source that is a force, not a moment
->> tensor. This is easy in the 2D code (elastic force option), but it
->> doesn't appear to be present in the 3D code. Can you confirm? Is that
->> something we want? (It seems like it.)
->>
->> Thanks,
->> Carl
-
--------- Original Message --------
-Subject: Re: force in SPECFEM3D
-Date: Thu, 19 Jul 2012 02:15:05 +0200
-From: Dimitri Komatitsch
-To: Carl Tape, daniel peter
-CC: Jeroen Tromp
-
-Hi Carl and Daniel,
-
-Yes, I think that would be very useful;
-and we should also put these parameters in the input file (Par_file)
-rather than in constants.h, so that users do not need to recompile the code.
-
-Thanks,
-Best wishes,
-
-Dimitri.
-
-On 07/19/2012 01:06 AM, Carl Tape wrote:
-> Hi Daniel,
->
-> Great. This is it:
->
-> integer, parameter :: COMPONENT_FORCE_SOURCE = 3 ! takes direction
-> in comp E/N/Z = 1/2/3
->
-> But perhaps we can generalize this to specify a unit vector that
-> doesn't have to be normalized, like
->
-> COMPONENT_FORCE_SOURCE_X = 1
-> COMPONENT_FORCE_SOURCE_Y = -2
-> COMPONENT_FORCE_SOURCE_Z = -1
->
-> Then in the code we could normalize (one less thing to worry about).
-> This would be closer to the 2D code (which takes and angle), and it
-> would allow for computing Green's functions on fault with arbitrary
-> orientation.
->
-> Let me know what you think.
->
-> Thanks,
-> Carl
->
-> On Wed, Jul 18, 2012 at 2:45 PM, daniel peter<dpeter at princeton.edu> wrote:
->> Hi Carl
->>
->> there is such an option in the constants.h file for the 3D codes (you will have to set the flag, something like "USE_POINT_FORCE" and direction, and then recompile the code).
->>
->> I'll see what we could put into the manual as well.
->>
->> Best wishes
->> Daniel
-
-------------------------------------------------
PT-SCOTCH:
------------------------------------------------
@@ -760,7 +526,240 @@
CPUs or GPUs, depending on the case; no mixed / hybrid model);
we then "svn delete" SUNFLOWER and keep a single code
+------------------------------------------------
+Par_file modifications:
+------------------------------------------------
+- suggestion 05:
+----------------
+
+Subject: Re: [CIG-SEISMO] simulation of seismic wave propagation excited by a force
+Date: Tue, 11 Sep 2012 18:34:21 +0200
+From: Dimitri Komatitsch
+To: Christina Morency
+CC: Yang Luo, Yingzi Ying
+
+Hi Christina and Yang, Hi all,
+
+OK, great. This is excellent news.
+
+So let us just move this "USE_FORCE_POINT_SOURCE" flag from
+setup/constants.h.in to DATA/Par_file, and let us mention it in the manual.
+
+Jo, could you please do that and commit the changes?
+
+Thank you,
+Dimitri.
+
+On 09/11/2012 06:17 PM, Christina Morency wrote:
+> Hi Dimitri,
+>
+> Yes, I can confirm what Yang is saying, it is already in the package,
+> with interpolation if the source is not located on a GLL point (as it is
+> indeed done in the 2D btw).
+>
+> Christina
+>
+> On 09/11/2012 09:14 AM, Yang Luo wrote:
+>> I think Daniel (or someone else) has already copied necessary lines from old package.
+>>
+>> There is an option "USE_FORCE_POINT_SOURCE" in constant.h, which allows one to use a point force
+>> (which Yingzi already figured out, or maybe it's already in the manual)
+>> logical, parameter :: USE_FORCE_POINT_SOURCE = .false.
+>> double precision, parameter :: FACTOR_FORCE_SOURCE = 1.d15
+>> integer, parameter :: COMPONENT_FORCE_SOURCE = 3 ! takes direction in comp E/N/Z = 1/2/3
+>>
+>> Interpolation is included, if the point force is not located exactly on a GLL point.
+
+=================
+
+add an option to use a point force source in SPECFEM3D instead of a CMTSOLUTION; some users need that (very classical for instance
+in non destructive testing) and it is easy to do; the only thing to decide is what convention to use in DATA/Par_file: add a flag
+for that?
+put a flag to choose between a Ricker source and a Heaviside source when a force source is used
+(we can cut and paste the Heaviside implementation from the case of a CMT source)
+
+Feedback from Qinya: I think this is a great idea. In the past, we always had multiple versions of the same specfem code, just to be able to run point force sources. We can add a parameter `SOURCE_TYPE' in Par file, which is = 0 for moment-tensor source ( reads CMTSOLUTION, and assumes Heavside/err as source time function), and = 1 for point force source (and read something like PNTSOLUTION, including force location, direction, time shift and half duration). Now there is no standard stf to use for point force, but Heavside function does not seem to be a good idea because it will produce a static offset on all seismograms, which is not desirable. On the other hand, I don't see why Richer wavelet (2nd order derivative of Gaussian) is necessary better than Gaussian itself. One way to give user the freedom of choice is to put the stf info also in the PNTSOLUTION.
+
+Subject: Re: [CIG-SEISMO] simulation of seismic wave propagation excited by a force
+Date: Tue, 11 Sep 2012 18:01:38 +0200
+From: Dimitri Komatitsch
+To: Yingzi Ying
+
+Dear Yingzi,
+
+You are right, in the current code there is no option in the Par_file to
+use a point source, but we added this to the todo list right before the
+summer and thus we should/will definitely do it. It is easy to do: just
+add the Ricker wavelet to the accel() vector (which is in fact a force
+vector before we divide it by the mass matrix), to the component you
+want (vertical or horizontal) and you are all set.
+
+Joseph, could you maybe do it? below are a few lines from one of my test
+codes, in which I do just that; let us cut and paste these lines in
+SPECFEM3D (it is just a matter of selecting the right grid point; which
+we already do for a CMT source, thus we could add a force instead of a
+CMT at that location; the only difference being that the CMT source is
+added to all the points of a given spectral element, while a force
+vector should be added to a single point).
+
+Let me talk to Joseph tomorrow to see if we can release a patch of the
+source code soon.
+
+I also cc Yang to see if he has already done that in a local version at
+Princeton maybe (for instance to model a force for some of his oil
+industry simulations in foothill regions?).
+
+Another difficulty is if the location of the force does not correspond
+to a grid point; then there are two options:
+
+- select the closest grid point; i.e. slightly change the location of
+the source (can be done automatically by the code)
+
+- use Christina Morency's way of adding a force source between grid
+points; I think I remember she did that in the 2D code, but I am not
+sure how; thus I cc her.
+
+(that second option would be better, if Christina can confirm that it
+works fine)
+
+Best wishes,
+Dimitri.
+
+==================================================================
+
+! ipoinsource is the global grid point at which we add the source
+! "3" means a vertical force (i.e. we add the force to the vertical
+! component of the force vector, which is later divided by the mass
+! matrix to get the acceleration vector)
+! factor_force is the norm of the force vector to add
+! ricker(t) is the Ricker wavelet at time t = (it-1)*deltat
+
+ accel(3,ipoinsource) = accel(3,ipoinsource) - factor_force * ricker(t)
+
+==================================================================
+
+On 09/11/2012 12:59 PM, Yingzi Ying wrote:
+> Hi all,
+>
+> I am trying to use specfem3d to simulate waves excited by a vertical
+> force with user-defined source forms such as ricker wavelet.
+>
+> I tried to do it through adjoint simulation by setting 'SIMULATION_TYPE
+> = 2' and put source position at STATIONS_ADJOINT and receivers at
+> CMTSOLUTION. But I failed.
+>
+> Does anyone have an idea or suggestion how to do active simulations with
+> specfem3d?
+>
+> Thanks in advance.
+>
+> Best regards,
+> Yingzi
+>
+
+
+-------- Original Message --------
+Subject: Re: force in SPECFEM3D
+Date: Wed, 18 Jul 2012 23:50:19 +0200
+From: Dimitri Komatitsch
+To: Carl Tape, daniel peter
+
+Hi all,
+
+Yes, let us definitely add it and make it an option in the Par_file.
+This will be very useful also in the case of non destructive testing of
+materials (in which case the source is almost always a set of forces).
+
+What probably exists is a routine to add a force at a GLL grid point
+(which is easy to do). If we also want to have an option to model a
+source located somewhere between grid points, then I guess we will need
+to see how the discrete source term looks like (I think Daniel once told
+me that Christina had done it in the 2D code; probably summing the
+Lagrange interpolants at all the GLL points of the element? (as we do
+for receivers located between grid points?).
+
+Thanks,
+Dimitri.
+
+On 07/18/2012 11:01 PM, Jeroen Tromp wrote:
+> Hi Daniel:
+>
+> The adjoint source involves a "force", so the necessary code exists. But
+> it can't hurt to explicitly allow a source that is a force rather than a
+> moment tensor.
+>
+> Best regards,
+>
+> Jeroen
+>
+> On 7/18/12 8:56 PM, Carl Tape wrote:
+>> Hi Daniel, Dimitri, Jeroen--
+>>
+>> I am continuing some benchmark simulations for static offsets in
+>> SPECFEM3D: 0-D, 1-D, 3-D models. One of the applications is to make a
+>> Green's function database similar to what is used in GPS inversions of
+>> fault slip. This means using a source that is a force, not a moment
+>> tensor. This is easy in the 2D code (elastic force option), but it
+>> doesn't appear to be present in the 3D code. Can you confirm? Is that
+>> something we want? (It seems like it.)
+>>
+>> Thanks,
+>> Carl
+
+-------- Original Message --------
+Subject: Re: force in SPECFEM3D
+Date: Thu, 19 Jul 2012 02:15:05 +0200
+From: Dimitri Komatitsch
+To: Carl Tape, daniel peter
+CC: Jeroen Tromp
+
+Hi Carl and Daniel,
+
+Yes, I think that would be very useful;
+and we should also put these parameters in the input file (Par_file)
+rather than in constants.h, so that users do not need to recompile the code.
+
+Thanks,
+Best wishes,
+
+Dimitri.
+
+On 07/19/2012 01:06 AM, Carl Tape wrote:
+> Hi Daniel,
+>
+> Great. This is it:
+>
+> integer, parameter :: COMPONENT_FORCE_SOURCE = 3 ! takes direction
+> in comp E/N/Z = 1/2/3
+>
+> But perhaps we can generalize this to specify a unit vector that
+> doesn't have to be normalized, like
+>
+> COMPONENT_FORCE_SOURCE_X = 1
+> COMPONENT_FORCE_SOURCE_Y = -2
+> COMPONENT_FORCE_SOURCE_Z = -1
+>
+> Then in the code we could normalize (one less thing to worry about).
+> This would be closer to the 2D code (which takes and angle), and it
+> would allow for computing Green's functions on fault with arbitrary
+> orientation.
+>
+> Let me know what you think.
+>
+> Thanks,
+> Carl
+>
+> On Wed, Jul 18, 2012 at 2:45 PM, daniel peter<dpeter at princeton.edu> wrote:
+>> Hi Carl
+>>
+>> there is such an option in the constants.h file for the 3D codes (you will have to set the flag, something like "USE_POINT_FORCE" and direction, and then recompile the code).
+>>
+>> I'll see what we could put into the manual as well.
+>>
+>> Best wishes
+>> Daniel
+
-----------------------------------------------------------
Stacey conditions:
-----------------------------------------------------------
More information about the CIG-COMMITS
mailing list