[cig-commits] r9275 - seismo/2D/SPECFEM2D/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Feb 10 15:37:12 PST 2008
Author: dkomati1
Date: 2008-02-10 15:37:12 -0800 (Sun, 10 Feb 2008)
New Revision: 9275
Modified:
seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
seismo/2D/SPECFEM2D/trunk/locate_source_force.F90
seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
Log:
cleaned the implementation and normalization of collocated force source
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2008-02-09 21:48:41 UTC (rev 9274)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2008-02-10 23:37:12 UTC (rev 9275)
@@ -564,32 +564,21 @@
endif ! end of absorbing boundaries
-! --- add the source
+! --- add the source if it is a moment tensor
if(.not. initialfield) then
! if this processor carries the source and the source element is elastic
if (is_proc_source == 1 .and. elastic(ispec_selected_source)) then
-! collocated force
-! beware, for acoustic medium, source is a potential, therefore source time function
-! gives shape of velocity, not displacement
- if(source_type == 1) then
-
- accel_elastic(1,iglob_source) = accel_elastic(1,iglob_source) - sin(angleforce)*source_time_function(it)
- accel_elastic(2,iglob_source) = accel_elastic(2,iglob_source) + cos(angleforce)*source_time_function(it)
-
! moment tensor
- else if(source_type == 2) then
+ if(source_type == 2) then
! add source array
- do j=1,NGLLZ
- do i=1,NGLLX
- iglob = ibool(i,j,ispec_selected_source)
- accel_elastic(:,iglob) = accel_elastic(:,iglob) + sourcearray(:,i,j)*source_time_function(it)
- enddo
- enddo
-
- else
- call exit_MPI('wrong source type in elastic element')
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source)
+ accel_elastic(:,iglob) = accel_elastic(:,iglob) + sourcearray(:,i,j)*source_time_function(it)
+ enddo
+ enddo
endif
endif ! if this processor carries the source and the source element is elastic
Modified: seismo/2D/SPECFEM2D/trunk/locate_source_force.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_source_force.F90 2008-02-09 21:48:41 UTC (rev 9274)
+++ seismo/2D/SPECFEM2D/trunk/locate_source_force.F90 2008-02-10 23:37:12 UTC (rev 9275)
@@ -40,8 +40,8 @@
!
!========================================================================
- subroutine locate_source_force(coord,ibool,npoin,nspec,x_source,z_source,source_type,ix_source,iz_source, &
- ispec_source,iglob_source,is_proc_source,nb_proc_source)
+ subroutine locate_source_force(coord,ibool,npoin,nspec,x_source,z_source,ix_source,iz_source, &
+ ispec_source,iglob_source,is_proc_source,nb_proc_source)
!
!----- calculer la position reelle de la source
@@ -54,7 +54,7 @@
include "mpif.h"
#endif
- integer npoin,nspec,source_type
+ integer npoin,nspec
integer ibool(NGLLX,NGLLZ,nspec)
double precision x_source,z_source
@@ -80,29 +80,22 @@
ihighx = NGLLX
ihighz = NGLLZ
-! on ne fait la recherche que sur l'interieur de l'element si source explosive
- if(source_type == 2) then
- ilowx = 2
- ilowz = 2
- ihighx = NGLLX-1
- ihighz = NGLLZ-1
- endif
+! look for the closest grid point
+ do numelem = 1,nspec
-! recherche du point de grille le plus proche
- do numelem=1,nspec
- do ix=ilowx,ihighx
- do iz=ilowz,ihighz
+ do ix = ilowx,ihighx
+ do iz = ilowz,ihighz
-! numero global du point
- ip=ibool(ix,iz,numelem)
+! global point number
+ ip = ibool(ix,iz,numelem)
-! coordonnees du point de grille
+! coordinates of this grid point
xp = coord(1,ip)
zp = coord(2,ip)
dist = sqrt((xp-x_source)**2 + (zp-z_source)**2)
-! retenir le point pour lequel l'ecart est minimal
+! keep the point for which distance is minimum
if(dist < distmin) then
distmin = dist
iglob_source = ip
@@ -113,13 +106,13 @@
enddo
enddo
+
enddo
distminmax = max(distmin,distminmax)
-
#ifdef USE_MPI
- ! global minimum distance computed over all processes
+! global minimum distance computed over all processes
call MPI_ALLREDUCE (distminmax, dist_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ierror)
#else
@@ -127,15 +120,11 @@
#endif
+! check if this process contains the source
+ if ( dist_glob == distminmax ) is_proc_source = 1
- ! check if this process contains the source
- if ( dist_glob == distminmax ) then
- is_proc_source = 1
- end if
-
-
#ifdef USE_MPI
- ! determining the number of processes that contain the source (useful when the source is located on an interface)
+! determining the number of processes that contain the source (useful when the source is located on an interface)
call MPI_ALLREDUCE (is_proc_source, nb_proc_source, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierror)
#else
Modified: seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90 2008-02-09 21:48:41 UTC (rev 9274)
+++ seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90 2008-02-10 23:37:12 UTC (rev 9275)
@@ -91,24 +91,24 @@
! **************
if ( myrank == 0 .or. nproc == 1 ) then
- write(IOUT,*)
- write(IOUT,*) '*******************************'
- write(IOUT,*) ' locating moment-tensor source'
- write(IOUT,*) '*******************************'
- write(IOUT,*)
+ write(IOUT,*)
+ write(IOUT,*) '*******************************'
+ write(IOUT,*) ' locating moment-tensor source'
+ write(IOUT,*) '*******************************'
+ write(IOUT,*)
end if
! set distance to huge initial value
- distmin=HUGEVAL
+ distmin = HUGEVAL
is_proc_source = 0
- do ispec=1,nspec
+ do ispec = 1,nspec
! loop only on points inside the element
! exclude edges to ensure this point is not shared with other elements
- do j=2,NGLLZ-1
- do i=2,NGLLX-1
+ do j = 2,NGLLZ-1
+ do i = 2,NGLLX-1
iglob = ibool(i,j,ispec)
dist = sqrt((x_source-dble(coord(1,iglob)))**2 + (z_source-dble(coord(2,iglob)))**2)
@@ -127,7 +127,6 @@
! end of loop on all the spectral elements
enddo
-
#ifdef USE_MPI
! global minimum distance computed over all processes
call MPI_ALLREDUCE (distmin, dist_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ierror)
@@ -137,13 +136,9 @@
#endif
+! check if this process contains the source
+ if ( dist_glob == distmin ) is_proc_source = 1
- ! check if this process contains the source
- if ( dist_glob == distmin ) then
- is_proc_source = 1
- end if
-
-
#ifdef USE_MPI
! determining the number of processes that contain the source (useful when the source is located on an interface)
call MPI_ALLREDUCE (is_proc_source, nb_proc_source, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierror)
@@ -175,8 +170,8 @@
! ****************************************
! use initial guess in xi and gamma
- xi = xigll(ix_initial_guess)
- gamma = zigll(iz_initial_guess)
+ xi = xigll(ix_initial_guess)
+ gamma = zigll(iz_initial_guess)
! iterate to solve the non linear system
do iter_loop = 1,NUM_ITER
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2008-02-09 21:48:41 UTC (rev 9274)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2008-02-10 23:37:12 UTC (rev 9275)
@@ -194,7 +194,7 @@
double precision, dimension(:), allocatable :: vp_display
double precision, dimension(:,:,:), allocatable :: vpext,vsext,rhoext
- double precision :: previous_vsext
+ double precision :: previous_vsext,rho_at_source_location
double precision, dimension(:,:,:), allocatable :: shape2D,shape2D_display
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: xix,xiz,gammax,gammaz,jacobian
@@ -495,6 +495,18 @@
endif
endif
+! if Dirac source time function, use a very thin Gaussian instead
+! if Heaviside source time function, use a very thin error function instead
+! time delay of the source in seconds, use a 20 % security margin (use 2 / f0 if error function)
+ if(time_function_type == 4 .or. time_function_type == 5) then
+ f0 = 1.d0 / (10.d0 * deltat)
+ if(time_function_type == 5) then
+ t0 = 2.0d0 / f0
+ else
+ t0 = 1.20d0 / f0
+ endif
+ endif
+
! for the source time function
aval = pi*pi*f0*f0
@@ -723,7 +735,6 @@
end if
-
!
!---- close input file
!
@@ -959,12 +970,21 @@
!---- define actual location of source and receivers
if(source_type == 1) then
+
! collocated force source
- call locate_source_force(coord,ibool,npoin,nspec,x_source,z_source,source_type, &
+ call locate_source_force(coord,ibool,npoin,nspec,x_source,z_source, &
ix_source,iz_source,ispec_selected_source,iglob_source,is_proc_source,nb_proc_source)
+! get density at the source in order to implement collocated force with the right
+! amplitude later
+ if(is_proc_source == 1) then
+ rho_at_source_location = density(kmato(ispec_selected_source))
+! external velocity model
+ if(assign_external_model) rho_at_source_location = rhoext(ix_source,iz_source,ispec_selected_source)
+ endif
+
! check that acoustic source is not exactly on the free surface because pressure is zero there
- if ( is_proc_source == 1 ) then
+ if(is_proc_source == 1) then
do ispec_acoustic_surface = 1,nelem_acoustic_surface
ispec = acoustic_surface(1,ispec_acoustic_surface)
if( .not. elastic(ispec) .and. ispec == ispec_selected_source ) then
@@ -993,7 +1013,6 @@
call exit_MPI('incorrect source type')
endif
-
! locate receivers in the mesh
call locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank,&
st_xval,st_zval,ispec_selected_rec, &
@@ -1220,7 +1239,6 @@
#endif
-
! fill mass matrix with fictitious non-zero values to make sure it can be inverted globally
if(any_elastic) where(rmass_inverse_elastic <= 0._CUSTOM_REAL) rmass_inverse_elastic = 1._CUSTOM_REAL
if(any_acoustic) where(rmass_inverse_acoustic <= 0._CUSTOM_REAL) rmass_inverse_acoustic = 1._CUSTOM_REAL
@@ -1237,8 +1255,6 @@
! convert receiver angle to radians
anglerec = anglerec * pi / 180.d0
-
-
!
!---- for color images
!
@@ -1670,10 +1686,10 @@
allocate(source_time_function(NSTEP))
if ( myrank == 0 ) then
- write(IOUT,*)
- write(IOUT,*) 'Saving the source time function in a text file...'
- write(IOUT,*)
- open(unit=55,file='OUTPUT_FILES/source.txt',status='unknown')
+ write(IOUT,*)
+ write(IOUT,*) 'Saving the source time function in a text file...'
+ write(IOUT,*)
+ open(unit=55,file='OUTPUT_FILES/source.txt',status='unknown')
endif
! loop on all the time steps
@@ -1705,14 +1721,10 @@
endif
! output absolute time in third column, in case user wants to check it as well
- if ( myrank == 0 ) then
- write(55,*) sngl(time),real(source_time_function(it),4),sngl(time-t0)
- endif
+ if (myrank == 0) write(55,*) sngl(time),real(source_time_function(it),4),sngl(time-t0)
enddo
- if ( myrank == 0 ) then
- close(55)
- endif
+ if (myrank == 0) close(55)
! nb_proc_source is the number of processes that own the source (the nearest point). It can be greater
! than one if the nearest point is on the interface between several partitions with an explosive source.
@@ -2313,13 +2325,33 @@
! ************************************************************************************
if(any_elastic) then
+
accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic
accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic
+
+! --- add the source if it is a collocated force
+ if(.not. initialfield) then
+
+! if this processor carries the source and the source element is elastic
+ if (is_proc_source == 1 .and. elastic(ispec_selected_source)) then
+
+! collocated force
+ if(source_type == 1) then
+ accel_elastic(1,iglob_source) = accel_elastic(1,iglob_source) &
+ - sin(angleforce)*source_time_function(it) / rho_at_source_location
+ accel_elastic(2,iglob_source) = accel_elastic(2,iglob_source) &
+ + cos(angleforce)*source_time_function(it) / rho_at_source_location
+ endif
+
+ endif ! if this processor carries the source and the source element is elastic
+
+ endif ! if not using an initial field
+
veloc_elastic = veloc_elastic + deltatover2*accel_elastic
+
endif
!---- compute kinetic and potential energy
-!
if(OUTPUT_ENERGY) &
call compute_energy(displ_elastic,veloc_elastic, &
xix,xiz,gammax,gammaz,jacobian,ibool,elastic,hprime_xx,hprime_zz, &
Modified: seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2008-02-09 21:48:41 UTC (rev 9274)
+++ seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2008-02-10 23:37:12 UTC (rev 9275)
@@ -58,7 +58,6 @@
integer :: NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current
double precision :: t0,deltat
-
integer, intent(in) :: nrecloc,myrank
integer, dimension(nrec),intent(in) :: which_proc_receiver
More information about the cig-commits
mailing list