[cig-commits] r13290 - seismo/3D/SPECFEM3D_SESAME/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Tue Nov 11 11:10:33 PST 2008
Author: dkomati1
Date: 2008-11-11 11:10:33 -0800 (Tue, 11 Nov 2008)
New Revision: 13290
Modified:
seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces.f90
seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in
seismo/3D/SPECFEM3D_SESAME/trunk/locate_source.f90
seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90
Log:
clarified the way a force source (e.g. vertical or normal) is entered
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces.f90 2008-11-11 17:21:09 UTC (rev 13289)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces.f90 2008-11-11 19:10:33 UTC (rev 13290)
@@ -263,7 +263,7 @@
if (ispec_is_inner(ispec_selected_source(isource)) .eqv. phase_is_inner) then
- if(FASTER_SOURCES_POINTS_ONLY) then
+ if(USE_FORCE_POINT_SOURCE) then
! add the source (only if this proc carries the source)
if(myrank == islice_selected_source(isource)) then
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in 2008-11-11 17:21:09 UTC (rev 13289)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in 2008-11-11 19:10:33 UTC (rev 13290)
@@ -124,8 +124,13 @@
! no lagrange interpolation on seismograms (we take the value on one NGLL point)
logical, parameter :: FASTER_RECEIVERS_POINTS_ONLY = .false.
- logical, parameter :: FASTER_SOURCES_POINTS_ONLY = .false.
+! 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.
+ logical, parameter :: USE_FORCE_POINT_SOURCE = .false.
+ double precision, parameter :: FACTOR_FORCE_SOURCE = 1.d10
+
! the receivers can be located inside the model
logical, parameter :: RECVS_CAN_BE_BURIED_EXT_MESH = .true.
logical, parameter :: SOURCES_CAN_BE_BURIED_EXT_MESH = .true.
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/locate_source.f90 2008-11-11 17:21:09 UTC (rev 13289)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/locate_source.f90 2008-11-11 19:10:33 UTC (rev 13290)
@@ -278,7 +278,7 @@
! define the interval in which we look for points
- if(FASTER_SOURCES_POINTS_ONLY .and. USE_EXTERNAL_MESH) then
+ if(USE_FORCE_POINT_SOURCE .and. USE_EXTERNAL_MESH) then
imin = 1
imax = NGLLX
@@ -490,7 +490,7 @@
! find the best (xi,eta,gamma) for the source
! *******************************************
- if(.not. FASTER_SOURCES_POINTS_ONLY) then
+ if(.not. USE_FORCE_POINT_SOURCE) then
! use initial guess in xi, eta and gamma
xi = xigll(ix_initial_guess_source)
@@ -587,7 +587,7 @@
final_distance_source(isource) = dsqrt((x_target_source-x_found_source(isource))**2 + &
(y_target_source-y_found_source(isource))**2 + (z_target_source-z_found_source(isource))**2)
- endif ! of if (.not. FASTER_SOURCES_POINTS_ONLY)
+ endif ! of if (.not. USE_FORCE_POINT_SOURCE)
! end of loop on all the sources
enddo
@@ -659,7 +659,7 @@
write(IMAIN,*) 'source located in slice ',islice_selected_source(isource)
write(IMAIN,*) ' in element ',ispec_selected_source(isource)
write(IMAIN,*)
- if(FASTER_SOURCES_POINTS_ONLY) then
+ if(USE_FORCE_POINT_SOURCE) then
write(IMAIN,*) ' xi coordinate of source in that element: ',nint(xi_source(isource))
write(IMAIN,*) ' eta coordinate of source in that element: ',nint(eta_source(isource))
write(IMAIN,*) 'gamma coordinate of source in that element: ',nint(gamma_source(isource))
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90 2008-11-11 17:21:09 UTC (rev 13289)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90 2008-11-11 19:10:33 UTC (rev 13290)
@@ -2890,11 +2890,13 @@
if (SIMULATION_TYPE == 1) then
-!!!!!!!!!! do isource = 1,NSOURCES
- do isource = 1,1
+! only one source for now if we use a force located at a grid point
+ if(USE_FORCE_POINT_SOURCE) NSOURCES = 1
- if(FASTER_SOURCES_POINTS_ONLY) then
+ do isource = 1,NSOURCES
+ if(USE_FORCE_POINT_SOURCE) then
+
! add the source (only if this proc carries the source)
if(myrank == islice_selected_source(isource)) then
@@ -2902,23 +2904,20 @@
nint(eta_source(isource)), &
nint(gamma_source(isource)), &
ispec_selected_source(isource))
+
f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
t0 = 1.2d0/f0
-if (it == 1 .and. myrank == 0) then
-
-print *,'using a source of dominant frequency ',f0
-print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
-print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
-endif
+
+ if (it == 1 .and. myrank == 0) then
+ print *,'using a source of dominant frequency ',f0
+ print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+ print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+ endif
-
- ! we use nu_source(:,3) here because we want a source normal to the surface.
- ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- !accel(:,iglob) = accel(:,iglob) + &
- ! sngl(nu_source(:,3,isource) * 10000000.d0 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
- ! exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
-accel(:,iglob) = accel(:,iglob) + &
- sngl(nu_source(:,3,isource) * 1.d10 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
+! we use nu_source(:,3) here because we want a source normal to the surface.
+! This is the expression of a Ricker
+ accel(:,iglob) = accel(:,iglob) + &
+ sngl(nu_source(:,3,isource) * FACTOR_FORCE_SOURCE * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
endif
@@ -2949,7 +2948,7 @@
endif
- endif ! end of if(FASTER_SOURCES_POINTS_ONLY)
+ endif ! end of if(USE_FORCE_POINT_SOURCE)
enddo
@@ -3051,11 +3050,13 @@
buffer_send_faces_vector,buffer_received_faces_vector,npoin2D_xi,npoin2D_eta, &
NPROC_XI,NPROC_ETA,NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY)
+! multiply by the inverse of the mass matrix
do i=1,NGLOB_AB
accel(1,i) = accel(1,i)*rmass(i)
accel(2,i) = accel(2,i)*rmass(i)
accel(3,i) = accel(3,i)*rmass(i)
enddo
+
if (SIMULATION_TYPE == 3) then
do i=1,NGLOB_AB
b_accel(1,i) = b_accel(1,i)*rmass(i)
More information about the CIG-COMMITS
mailing list