[cig-commits] r22507 - seismo/3D/SPECFEM3D/trunk/src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Thu Jul 4 03:28:15 PDT 2013
Author: danielpeter
Date: 2013-07-04 03:28:14 -0700 (Thu, 04 Jul 2013)
New Revision: 22507
Modified:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic_calling_routine.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
Log:
bug fix for point forces (FORCESOLUTION), uses now exact positions also for point forces (behaviour can be changed by setting flag USE_BEST_LOCATION to .false. in locate_source.f90
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2013-07-03 15:55:46 UTC (rev 22506)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2013-07-04 10:28:14 UTC (rev 22507)
@@ -29,7 +29,6 @@
subroutine compute_add_sources_acoustic(NSPEC_AB,NGLOB_AB,potential_dot_dot_acoustic, &
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
- xi_source,eta_source,gamma_source, &
hdur,hdur_gaussian,tshift_src,dt,t0, &
sourcearrays,kappastore,ispec_is_acoustic,&
SIMULATION_TYPE,NSTEP, &
@@ -40,7 +39,7 @@
xigll,yigll,zigll,xi_receiver,eta_receiver,gamma_receiver,&
station_name,network_name,adj_source_file,nrec_local,number_receiver_global, &
pm1_source_encoding,nsources_local,USE_FORCE_POINT_SOURCE, &
- USE_RICKER_TIME_FUNCTION,factor_force_source
+ USE_RICKER_TIME_FUNCTION
implicit none
include "constants.h"
@@ -62,7 +61,6 @@
! source
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(NSOURCES) :: hdur,hdur_gaussian,hdur_tiny,tshift_src
double precision :: dt,t0
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrays
@@ -126,12 +124,6 @@
if(USE_FORCE_POINT_SOURCE) then
- ! note: for use_force_point_source xi/eta/gamma are in the range [1,NGLL*]
- iglob = ibool(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec)
-
f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing FORCESOLUTION file format
!if (it == 1 .and. myrank == 0) then
@@ -141,11 +133,9 @@
!endif
if( USE_RICKER_TIME_FUNCTION ) then
- stf_used = factor_force_source(isource) * sourcearrays(isource,1,i,j,k) * &
- comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_src(isource),f0)
+ stf_used = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_src(isource),f0)
else
- stf_used = factor_force_source(isource) * sourcearrays(isource,1,i,j,k) * &
- comp_source_time_function_gauss(dble(it-1)*DT-t0-tshift_src(isource),hdur_tiny(isource))
+ stf_used = comp_source_time_function_gauss(dble(it-1)*DT-t0-tshift_src(isource),hdur_tiny(isource))
endif
! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
@@ -154,10 +144,15 @@
! acoustic source for pressure gets divided by kappa
! source contribution
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - stf_used / kappastore(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)),ispec)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - sourcearrays(isource,1,i,j,k) * stf_used / kappastore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
else
@@ -373,7 +368,6 @@
subroutine compute_add_sources_acoustic_bpwf(NSPEC_AB, &
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
- xi_source,eta_source,gamma_source, &
hdur,hdur_gaussian,tshift_src,dt,t0, &
sourcearrays,kappastore,ispec_is_acoustic,&
SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
@@ -381,7 +375,7 @@
use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total, &
pm1_source_encoding,nsources_local,USE_FORCE_POINT_SOURCE, &
- USE_RICKER_TIME_FUNCTION,factor_force_source
+ USE_RICKER_TIME_FUNCTION
implicit none
include "constants.h"
@@ -400,7 +394,6 @@
! source
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(NSOURCES) :: hdur,hdur_gaussian,hdur_tiny,tshift_src
double precision :: dt,t0
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrays
@@ -474,14 +467,8 @@
if(USE_FORCE_POINT_SOURCE) then
- ! note: for use_force_point_source xi/eta/gamma are in the range [1,NGLL*]
- iglob = ibool(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec)
+ f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing FORCESOLUTION file format
- f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
-
!if (it == 1 .and. myrank == 0) then
! write(IMAIN,*) 'using a source of dominant frequency ',f0
! write(IMAIN,*) 'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
@@ -489,23 +476,28 @@
!endif
if( USE_RICKER_TIME_FUNCTION ) then
- stf_used = factor_force_source(isource) * sourcearrays(isource,1,i,j,k) * &
- comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_src(isource),f0)
+ stf_used = comp_source_time_function_rickr( &
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),f0)
else
- stf_used = factor_force_source(isource) * sourcearrays(isource,1,i,j,k) * &
- comp_source_time_function_gauss(dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur_tiny(isource))
+ stf_used = comp_source_time_function_gauss( &
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur_tiny(isource))
endif
! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
! the sign is negative because pressure p = - Chi_dot_dot therefore we need
! to add minus the source to Chi_dot_dot to get plus the source in pressure:
-
+
! acoustic source for pressure gets divided by kappa
! source contribution
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
- - stf_used / kappastore(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)),ispec)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
+ - sourcearrays(isource,1,i,j,k) * stf_used / kappastore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
else
if( USE_RICKER_TIME_FUNCTION ) then
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2013-07-03 15:55:46 UTC (rev 22506)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2013-07-04 10:28:14 UTC (rev 22507)
@@ -32,7 +32,6 @@
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, &
hdur,hdur_gaussian,tshift_src,dt,t0,sourcearrays, &
ispec_is_poroelastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
@@ -66,7 +65,6 @@
! source
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(NSOURCES) :: hdur,hdur_gaussian,tshift_src
double precision :: dt,t0
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrays
@@ -118,31 +116,6 @@
if(USE_FORCE_POINT_SOURCE) then
- ! note: for use_force_point_source xi/eta/gamma are in the range [1,NGLL*]
- iglob = ibool(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec_selected_source(isource))
-
- ! get poroelastic parameters of current local GLL
- phil = phistore(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec_selected_source(isource))
- tortl = tortstore(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec_selected_source(isource))
- rhol_s = rhoarraystore(1,nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec_selected_source(isource))
- rhol_f = rhoarraystore(2,nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec_selected_source(isource))
- rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
-
!f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
!t0 = 1.2d0/f0
@@ -170,17 +143,23 @@
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) + &
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ ! get poroelastic parameters of current local GLL
+ phil = phistore(i,j,k,ispec)
+ tortl = tortstore(i,j,k,ispec)
+ rhol_s = rhoarraystore(1,i,j,k,ispec)
+ rhol_f = rhoarraystore(2,i,j,k,ispec)
+ rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
+ ! solid phase
+ accels(:,iglob) = accels(:,iglob) + &
(1._CUSTOM_REAL - phil/tortl) * sourcearrays(isource,:,i,j,k) * stf_used
-! fluid phase
- accelw(:,iglob) = accelw(:,iglob) + &
+ ! fluid phase
+ accelw(:,iglob) = accelw(:,iglob) + &
(1._CUSTOM_REAL - rhol_f/rhol_bar) * sourcearrays(isource,:,i,j,k) * stf_used
- enddo
- enddo
+ enddo
+ enddo
enddo
else
@@ -190,36 +169,36 @@
!stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_src(isource),hdur(isource))
stf = comp_source_time_function_gauss(dble(it-1)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
- ! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- stf_used = sngl(stf)
- else
- stf_used = stf
- endif
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
- ! add source array
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
- ! get poroelastic parameters of current local GLL
- phil = phistore(i,j,k,ispec)
- tortl = tortstore(i,j,k,ispec)
- rhol_s = rhoarraystore(1,i,j,k,ispec)
- rhol_f = rhoarraystore(2,i,j,k,ispec)
- rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
-! source in the solid phase only
-! solid phase
- accels(:,iglob) = accels(:,iglob) &
+ ! add source array
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ ! get poroelastic parameters of current local GLL
+ phil = phistore(i,j,k,ispec)
+ tortl = tortstore(i,j,k,ispec)
+ rhol_s = rhoarraystore(1,i,j,k,ispec)
+ rhol_f = rhoarraystore(2,i,j,k,ispec)
+ rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
+ ! source in the solid phase only
+ ! solid phase
+ accels(:,iglob) = accels(:,iglob) &
+ sourcearrays(isource,:,i,j,k)*stf_used
! + (1._CUSTOM_REAL - phil/tortl)*sourcearrays(isource,:,i,j,k)*stf_used
-! fluid phase
- accelw(:,iglob) = accelw(:,iglob) &
+ ! fluid phase
+ accelw(:,iglob) = accelw(:,iglob) &
- rhol_f/rhol_bar*sourcearrays(isource,:,i,j,k)*stf_used
! + (1._CUSTOM_REAL - rhol_f/rhol_bar)*sourcearrays(isource,:,i,j,k)*stf_used
- enddo
enddo
- enddo
+ enddo
+ enddo
endif ! USE_FORCE_POINT_SOURCE
@@ -382,31 +361,6 @@
if(USE_FORCE_POINT_SOURCE) then
- ! note: for use_force_point_source xi/eta/gamma are in the range [1,NGLL*]
- iglob = ibool(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec_selected_source(isource))
-
- ! get poroelastic parameters of current local GLL
- phil = phistore(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec_selected_source(isource))
- tortl = tortstore(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec_selected_source(isource))
- rhol_s = rhoarraystore(1,nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec_selected_source(isource))
- rhol_f = rhoarraystore(2,nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec_selected_source(isource))
- rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
-
!f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
!if (it == 1 .and. myrank == 0) then
! write(IMAIN,*) 'using a source of dominant frequency ',f0
@@ -432,17 +386,23 @@
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) + &
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ ! get poroelastic parameters of current local GLL
+ phil = phistore(i,j,k,ispec)
+ tortl = tortstore(i,j,k,ispec)
+ rhol_s = rhoarraystore(1,i,j,k,ispec)
+ rhol_f = rhoarraystore(2,i,j,k,ispec)
+ rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
+ ! solid phase
+ 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) + &
+ ! fluid phase
+ b_accelw(:,iglob) = b_accelw(:,iglob) + &
(1._CUSTOM_REAL - rhol_f/rhol_bar) * sourcearrays(isource,:,i,j,k) * stf_used
- enddo
- enddo
+ enddo
+ enddo
enddo
else
@@ -465,20 +425,20 @@
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
- iglob = ibool(i,j,k,ispec_selected_source(isource))
- ! get poroelastic parameters of current local GLL
- phil = phistore(i,j,k,ispec_selected_source(isource))
- tortl = tortstore(i,j,k,ispec_selected_source(isource))
- rhol_s = rhoarraystore(1,i,j,k,ispec_selected_source(isource))
- rhol_f = rhoarraystore(2,i,j,k,ispec_selected_source(isource))
- rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
-! source in the solid phase only
-! solid phase
- b_accels(:,iglob) = b_accels(:,iglob) &
+ iglob = ibool(i,j,k,ispec)
+ ! get poroelastic parameters of current local GLL
+ phil = phistore(i,j,k,ispec)
+ tortl = tortstore(i,j,k,ispec)
+ rhol_s = rhoarraystore(1,i,j,k,ispec)
+ rhol_f = rhoarraystore(2,i,j,k,ispec)
+ rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
+ ! source in the solid phase only
+ ! solid phase
+ b_accels(:,iglob) = b_accels(:,iglob) &
+ sourcearrays(isource,:,i,j,k)*stf_used
! + (1._CUSTOM_REAL - phil/tortl)*sourcearrays(isource,:,i,j,k)*stf_used
-! fluid phase
- b_accelw(:,iglob) = b_accelw(:,iglob) &
+ ! fluid phase
+ b_accelw(:,iglob) = b_accelw(:,iglob) &
- rhol_f/rhol_bar*sourcearrays(isource,:,i,j,k)*stf_used
! + (1._CUSTOM_REAL - rhol_f/rhol_bar)*sourcearrays(isource,:,i,j,k)*stf_used
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-07-03 15:55:46 UTC (rev 22506)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-07-04 10:28:14 UTC (rev 22507)
@@ -177,8 +177,7 @@
call compute_add_sources_acoustic(NSPEC_AB,NGLOB_AB,potential_dot_dot_acoustic, &
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
- xi_source,eta_source,gamma_source, &
- hdur,hdur_gaussian,tshift_src,dt,t0, &
+ hdur,hdur_gaussian,tshift_src,dt,t0, &
sourcearrays,kappastore,ispec_is_acoustic,&
SIMULATION_TYPE,NSTEP, &
nrec,islice_selected_rec,ispec_selected_rec, &
@@ -421,7 +420,6 @@
call compute_add_sources_acoustic_bpwf(NSPEC_AB, &
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
- xi_source,eta_source,gamma_source, &
hdur,hdur_gaussian,tshift_src,dt,t0, &
sourcearrays,kappastore,ispec_is_acoustic,&
SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic_calling_routine.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic_calling_routine.f90 2013-07-03 15:55:46 UTC (rev 22506)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic_calling_routine.f90 2013-07-04 10:28:14 UTC (rev 22507)
@@ -205,7 +205,6 @@
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, &
hdur,hdur_gaussian,tshift_src,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/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2013-07-03 15:55:46 UTC (rev 22506)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2013-07-04 10:28:14 UTC (rev 22507)
@@ -161,6 +161,15 @@
integer, dimension(NSOURCES) :: idomain
integer, dimension(NGATHER_SOURCES,0:NPROC-1) :: idomain_all
+ !-----------------------------------------------------------------------------------
+
+ ! interpolates position to find best possible location closest to target position (default)
+ ! note: this can be turned off to speedup location search
+ ! and position source at closest GLL point
+ logical, parameter :: USE_BEST_LOCATION = .true.
+
+ !-----------------------------------------------------------------------------------
+
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', OUTPUT_FILES_PATH(1:len_trim(OUTPUT_FILES_PATH)))
@@ -267,9 +276,9 @@
distmin = HUGEVAL
ispec_selected_source(isource) = 0
- ix_initial_guess_source = 0
- iy_initial_guess_source = 0
- iz_initial_guess_source = 0
+ ix_initial_guess_source = 1
+ iy_initial_guess_source = 1
+ iz_initial_guess_source = 1
do ispec=1,NSPEC_AB
! define the interval in which we look for points
@@ -320,8 +329,7 @@
iz_initial_guess_source = k
! store xi,eta,gamma and x,y,z of point found
- ! note: they have range [1.0d0,NGLLX/Y/Z], used for point sources
- ! see e.g. in compute_add_source_elastic.f90
+ ! note: here they have range [1.0d0,NGLLX/Y/Z]
xi_source(isource) = dble(ix_initial_guess_source)
eta_source(isource) = dble(iy_initial_guess_source)
gamma_source(isource) = dble(iz_initial_guess_source)
@@ -504,9 +512,8 @@
! find the best (xi,eta,gamma) for the source
! *******************************************
- ! for point sources, the location will be exactly at a GLL point
- ! otherwise this tries to find best location
- if(.not. USE_FORCE_POINT_SOURCE) then
+ ! this tries to find best location
+ if( USE_BEST_LOCATION ) then
! uses actual location interpolators, in range [-1,1]
xi = xigll(ix_initial_guess_source)
@@ -597,6 +604,7 @@
xi_source(isource) = xi
eta_source(isource) = eta
gamma_source(isource) = gamma
+
x_found_source(isource) = x
y_found_source(isource) = y
z_found_source(isource) = z
@@ -604,9 +612,16 @@
! compute final distance between asked and found (converted to km)
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)
+ else
- endif ! of if (.not. USE_FORCE_POINT_SOURCE)
+ ! takes initial GLL point guess and uses actual location interpolators, in range [-1,1]
+ ! note: xi/eta/gamma will be in range [-1,1]
+ xi_source(isource) = xigll(ix_initial_guess_source)
+ eta_source(isource) = yigll(iy_initial_guess_source)
+ gamma_source(isource) = zigll(iz_initial_guess_source)
+ endif ! USE_BEST_LOCATION
+
! end of loop on all the sources
enddo
@@ -733,9 +748,11 @@
write(IMAIN,*)
if(USE_FORCE_POINT_SOURCE) then
- write(IMAIN,*) ' i index of source in that element: ',nint(xi_source(isource))
- 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,*) 'using force point source: '
+ write(IMAIN,*) ' xi coordinate of source in that element: ',xi_source(isource)
+ write(IMAIN,*) ' eta coordinate of source in that element: ',eta_source(isource)
+ write(IMAIN,*) ' gamma coordinate of source in that element: ',gamma_source(isource)
+
write(IMAIN,*)
write(IMAIN,*) ' component of direction vector in East direction: ',comp_dir_vect_source_E(isource)
write(IMAIN,*) ' component of direction vector in North direction: ',comp_dir_vect_source_N(isource)
@@ -760,6 +777,7 @@
write(IMAIN,*)
write(IMAIN,*) ' half duration -> frequency: ',hdur(isource),' seconds**(-1)'
else
+ write(IMAIN,*) 'using moment tensor source: '
write(IMAIN,*) ' xi coordinate of source in that element: ',xi_source(isource)
write(IMAIN,*) ' eta coordinate of source in that element: ',eta_source(isource)
write(IMAIN,*) ' gamma coordinate of source in that element: ',gamma_source(isource)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2013-07-03 15:55:46 UTC (rev 22506)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2013-07-04 10:28:14 UTC (rev 22507)
@@ -284,47 +284,24 @@
izmin = minval( free_surface_ijk(3,:,iface) )
izmax = maxval( free_surface_ijk(3,:,iface) )
- if( .not. USE_FORCE_POINT_SOURCE ) then
- ! xmin face
- if(ixmin==1 .and. ixmax==1) then
- if( xi_source(isource) < -0.99d0) is_on = .true.
- ! xmax face
- else if(ixmin==NGLLX .and. ixmax==NGLLX) then
- if( xi_source(isource) > 0.99d0) is_on = .true.
- ! ymin face
- else if(iymin==1 .and. iymax==1) then
- if( eta_source(isource) < -0.99d0) is_on = .true.
- ! ymax face
- else if(iymin==NGLLY .and. iymax==NGLLY) then
- if( eta_source(isource) > 0.99d0) is_on = .true.
- ! zmin face
- else if(izmin==1 .and. izmax==1 ) then
- if( gamma_source(isource) < -0.99d0) is_on = .true.
- ! zmax face
- else if(izmin==NGLLZ .and. izmax==NGLLZ ) then
- if( gamma_source(isource) > 0.99d0) is_on = .true.
- endif
- else
- ! note: for use_force_point_source xi/eta/gamma_source values are in the range [1,NGLL*]
- ! xmin face
- if(ixmin==1 .and. ixmax==1) then
- if( nint(xi_source(isource)) == 1) is_on = .true.
- ! xmax face
- else if(ixmin==NGLLX .and. ixmax==NGLLX) then
- if( nint(xi_source(isource)) == NGLLX) is_on = .true.
- ! ymin face
- else if(iymin==1 .and. iymax==1) then
- if( nint(eta_source(isource)) == 1) is_on = .true.
- ! ymax face
- else if(iymin==NGLLY .and. iymax==NGLLY) then
- if( nint(eta_source(isource)) == NGLLY) is_on = .true.
- ! zmin face
- else if(izmin==1 .and. izmax==1 ) then
- if( nint(gamma_source(isource)) == 1) is_on = .true.
- ! zmax face
- else if(izmin==NGLLZ .and. izmax==NGLLZ ) then
- if( nint(gamma_source(isource)) ==NGLLZ) is_on = .true.
- endif
+ ! xmin face
+ if(ixmin==1 .and. ixmax==1) then
+ if( xi_source(isource) < -0.99d0) is_on = .true.
+ ! xmax face
+ else if(ixmin==NGLLX .and. ixmax==NGLLX) then
+ if( xi_source(isource) > 0.99d0) is_on = .true.
+ ! ymin face
+ else if(iymin==1 .and. iymax==1) then
+ if( eta_source(isource) < -0.99d0) is_on = .true.
+ ! ymax face
+ else if(iymin==NGLLY .and. iymax==NGLLY) then
+ if( eta_source(isource) > 0.99d0) is_on = .true.
+ ! zmin face
+ else if(izmin==1 .and. izmax==1 ) then
+ if( gamma_source(isource) < -0.99d0) is_on = .true.
+ ! zmax face
+ else if(izmin==NGLLZ .and. izmax==NGLLZ ) then
+ if( gamma_source(isource) > 0.99d0) is_on = .true.
endif
endif ! free_surface_ispec
@@ -543,26 +520,22 @@
integer :: irec
integer :: i,j,k
integer :: icomp,itime,nadj_files_found,nadj_files_found_tot,ier
+! integer :: iglob
character(len=3),dimension(NDIM) :: comp
character(len=256) :: filename
- double precision :: hlagrange
-
- double precision, dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrayd
double precision, dimension(NGLLX) :: hxis,hpxis
double precision, dimension(NGLLY) :: hetas,hpetas
double precision, dimension(NGLLZ) :: hgammas,hpgammas
- double precision, dimension(:,:), allocatable :: hxis_store,hetas_store,hgammas_store
+ double precision, dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrayd
+ double precision :: hlagrange
! forward simulations
if (SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3) then
allocate(sourcearray(NDIM,NGLLX,NGLLY,NGLLZ), &
- sourcearrays(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ), &
- hxis_store(NSOURCES,NGLLX), &
- hetas_store(NSOURCES,NGLLY), &
- hgammas_store(NSOURCES,NGLLZ), stat=ier)
+ sourcearrays(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
if( ier /= 0 ) stop 'error allocating array sourcearray'
! compute source arrays
@@ -582,51 +555,53 @@
call lagrange_any(eta_source(isource),NGLLY,yigll,hetas,hpetas)
call lagrange_any(gamma_source(isource),NGLLZ,zigll,hgammas,hpgammas)
- hxis_store(isource,:) = hxis(:)
- hetas_store(isource,:) = hetas(:)
- hgammas_store(isource,:) = hgammas(:)
-
if (USE_FORCE_POINT_SOURCE) then ! use of FORCESOLUTION files
- ! note: for use_force_point_source xi/eta/gamma are in the range [1,NGLL*]
+ ! note: for use_force_point_source xi/eta/gamma are also in the range [-1,1], for exact positioning
- ! initializes source array
- sourcearray(:,:,:,:) = 0._CUSTOM_REAL
- sourcearrayd(:,:,:,:) = 0.d0
+ ! initializes source array
+ sourcearrayd(:,:,:,:) = 0.0d0
- ! calculates source array for interpolated location
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- hlagrange = hxis_store(isource,i) * hetas_store(isource,j) * hgammas_store(isource,k)
+ ! calculates source array for interpolated location
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ hlagrange = hxis(i) * hetas(j) * hgammas(k)
- ! acoustic source
- ! identical source array components in x,y,z-direction
- if( ispec_is_acoustic(ispec) ) then
- sourcearrayd(:,i,j,k) = hlagrange
- endif
+ ! acoustic source
+ ! identical source array components in x,y,z-direction
+ if( ispec_is_acoustic(ispec) ) then
+ sourcearrayd(:,i,j,k) = factor_force_source(isource) * hlagrange
+ endif
- ! elastic source
- if( ispec_is_elastic(ispec) ) then
- ! 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:
- sourcearrayd(:,i,j,k) = factor_force_source(isource) * hlagrange * &
- ( nu_source(1,:,isource) * comp_dir_vect_source_E(isource) + &
- nu_source(2,:,isource) * comp_dir_vect_source_N(isource) + &
- nu_source(3,:,isource) * comp_dir_vect_source_Z_UP(isource) ) / &
- sqrt( comp_dir_vect_source_E(isource)**2 + comp_dir_vect_source_N(isource)**2 + &
- comp_dir_vect_source_Z_UP(isource)**2 )
- endif
- enddo
+ ! elastic source
+ if( ispec_is_elastic(ispec) ) then
+ ! checks norm of component vector
+ if( sqrt( comp_dir_vect_source_E(isource)**2 &
+ + comp_dir_vect_source_N(isource)**2 &
+ + comp_dir_vect_source_Z_UP(isource)**2 ) < TINYVAL ) then
+ call exit_MPI(myrank,'error force point source: component vector has (almost) zero norm')
+ endif
+
+ ! 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:
+ sourcearrayd(:,i,j,k) = factor_force_source(isource) * hlagrange * &
+ ( nu_source(1,:,isource) * comp_dir_vect_source_E(isource) + &
+ nu_source(2,:,isource) * comp_dir_vect_source_N(isource) + &
+ nu_source(3,:,isource) * comp_dir_vect_source_Z_UP(isource) ) / &
+ sqrt( comp_dir_vect_source_E(isource)**2 + comp_dir_vect_source_N(isource)**2 + &
+ comp_dir_vect_source_Z_UP(isource)**2 )
+ endif
enddo
- enddo
+ enddo
+ enddo
- ! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- sourcearray(:,:,:,:) = sngl(sourcearrayd(:,:,:,:))
- else
- sourcearray(:,:,:,:) = sourcearrayd(:,:,:,:)
- endif
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ sourcearray(:,:,:,:) = sngl(sourcearrayd(:,:,:,:))
+ else
+ sourcearray(:,:,:,:) = sourcearrayd(:,:,:,:)
+ endif
else ! use of CMTSOLUTION files
@@ -661,6 +636,7 @@
endif
endif
+
! stores source excitations
sourcearrays(isource,:,:,:,:) = sourcearray(:,:,:,:)
More information about the CIG-COMMITS
mailing list