[cig-commits] r8445 - seismo/2D/SPECFEM2D/trunk
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:46:49 PST 2007
Author: walter
Date: 2007-12-07 15:46:47 -0800 (Fri, 07 Dec 2007)
New Revision: 8445
Added:
seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90
seismo/2D/SPECFEM2D/trunk/locate_source_force.f90
seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.f90
Removed:
seismo/2D/SPECFEM2D/trunk/positsource.f90
Modified:
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90
seismo/2D/SPECFEM2D/trunk/defarrays.f90
seismo/2D/SPECFEM2D/trunk/locate_receivers.f90
seismo/2D/SPECFEM2D/trunk/meshfem2D.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.f90
seismo/2D/SPECFEM2D/trunk/write_seismograms.f90
Log:
added general moment-tensor source to 2D code
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-07 23:46:47 UTC (rev 8444)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-07 23:46:47 UTC (rev 8445)
@@ -13,8 +13,8 @@
# Intel Linux
F90 = ifort
-#FLAGS=-O0 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check bounds
-FLAGS=-O3 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check nobounds
+FLAGS=-O0 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check bounds
+#FLAGS=-O3 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check nobounds
#
# g95 (free f95 compiler from http://www.g95.org, still under development, but works)
@@ -33,9 +33,10 @@
OBJS_SPECFEM2D = $O/checkgrid.o $O/datim.o $O/defarrays.o\
$O/lagrange_poly.o $O/gmat01.o $O/gll_library.o $O/plotgll.o $O/define_derivative_matrices.o\
- $O/plotpost.o $O/locate_receivers.o $O/positsource.o $O/compute_gradient_attenuation.o\
+ $O/plotpost.o $O/locate_receivers.o $O/locate_source_force.o $O/compute_gradient_attenuation.o\
$O/specfem2D.o $O/write_seismograms.o $O/createnum_fast.o $O/createnum_slow.o\
- $O/define_shape_functions.o $O/cree_image_PNM.o $O/compute_gradient_fluid.o $O/recompute_jacobian.o
+ $O/define_shape_functions.o $O/cree_image_PNM.o $O/compute_gradient_fluid.o\
+ $O/recompute_jacobian.o $O/compute_arrays_source.o $O/locate_source_moment_tensor.o
default: meshfem2D specfem2D
@@ -98,9 +99,12 @@
$O/recompute_jacobian.o: recompute_jacobian.f90 constants.h
${F90} $(FLAGS) -c -o $O/recompute_jacobian.o recompute_jacobian.f90
-$O/positsource.o: positsource.f90 constants.h
- ${F90} $(FLAGS) -c -o $O/positsource.o positsource.f90
+$O/locate_source_force.o: locate_source_force.f90 constants.h
+ ${F90} $(FLAGS) -c -o $O/locate_source_force.o locate_source_force.f90
+$O/locate_source_moment_tensor.o: locate_source_moment_tensor.f90 constants.h
+ ${F90} $(FLAGS) -c -o $O/locate_source_moment_tensor.o locate_source_moment_tensor.f90
+
$O/define_shape_functions.o: define_shape_functions.f90 constants.h
${F90} $(FLAGS) -c -o $O/define_shape_functions.o define_shape_functions.f90
@@ -113,6 +117,9 @@
$O/compute_gradient_fluid.o: compute_gradient_fluid.f90 constants.h
${F90} $(FLAGS) -c -o $O/compute_gradient_fluid.o compute_gradient_fluid.f90
+$O/compute_arrays_source.o: compute_arrays_source.f90 constants.h
+ ${F90} $(FLAGS) -c -o $O/compute_arrays_source.o compute_arrays_source.f90
+
$O/cree_image_PNM.o: cree_image_PNM.f90 constants.h
${F90} $(FLAGS) -c -o $O/cree_image_PNM.o cree_image_PNM.f90
Added: seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90 2007-12-07 23:46:47 UTC (rev 8444)
+++ seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90 2007-12-07 23:46:47 UTC (rev 8445)
@@ -0,0 +1,99 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.1
+! ------------------------------
+!
+! Dimitri Komatitsch
+! Universite de Pau et des Pays de l'Adour, France
+!
+! (c) January 2005
+!
+!========================================================================
+
+ subroutine compute_arrays_source(ispec_selected_source,xi_source,gamma_source,sourcearray, &
+ Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,nspec)
+
+ implicit none
+
+ include "constants.h"
+
+ integer ispec_selected_source
+ integer nspec
+
+ double precision xi_source,gamma_source
+!!!!!!!!!!!!!!!!! double precision Mxx,Myy,Mzz,Mxy,Mxz,Myz
+ double precision Mxx,Mzz,Mxz
+
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
+
+ double precision, dimension(NDIM,NGLLX,NGLLZ) :: sourcearray
+
+ double precision xixd,xizd,gammaxd,gammazd
+
+! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: xigll
+ double precision, dimension(NGLLZ) :: zigll
+
+! source arrays
+ double precision, dimension(NGLLX,NGLLZ) :: G11,G13,G31,G33
+ double precision, dimension(NGLLX) :: hxis,hpxis
+ double precision, dimension(NGLLZ) :: hgammas,hpgammas
+
+ integer k,m
+ integer ir,iv
+
+! calculate G_ij for general source location
+! the source does not necessarily correspond to a Gauss-Lobatto point
+ do m=1,NGLLZ
+ do k=1,NGLLX
+
+ xixd = xix(k,m,ispec_selected_source)
+ xizd = xiz(k,m,ispec_selected_source)
+ gammaxd = gammax(k,m,ispec_selected_source)
+ gammazd = gammaz(k,m,ispec_selected_source)
+
+ G11(k,m) = Mxx*xixd+Mxz*xizd
+ G13(k,m) = Mxx*gammaxd+Mxz*gammazd
+ G31(k,m) = Mxz*xixd+Mzz*xizd
+ G33(k,m) = Mxz*gammaxd+Mzz*gammazd
+
+!!!! G21(k,m) = Mxy*xixd+Myz*xizd
+!!!! G23(k,m) = Mxy*gammaxd+Myz*gammazd
+
+ enddo
+ enddo
+
+! compute Lagrange polynomials at the source location
+ call lagrange_any(xi_source,NGLLX,xigll,hxis,hpxis)
+ call lagrange_any(gamma_source,NGLLZ,zigll,hgammas,hpgammas)
+
+! calculate source array
+ do m=1,NGLLZ
+ do k=1,NGLLX
+
+ sourcearray(:,k,m) = ZERO
+
+ do iv=1,NGLLZ
+ do ir=1,NGLLX
+
+ sourcearray(1,k,m) = sourcearray(1,k,m) + hxis(ir)*hgammas(iv) &
+ *(G11(ir,iv)*hpxis(k)*hgammas(m) &
+ +G13(ir,iv)*hxis(k)*hpgammas(m))
+
+! sourcearray(2,k,m) = sourcearray(2,k,m) + hxis(ir)*hgammas(iv) &
+! *(G21(ir,iv)*hpxis(k)*hgammas(m) &
+! +G23(ir,iv)*hxis(k)*hpgammas(m))
+
+ sourcearray(2,k,m) = sourcearray(2,k,m) + hxis(ir)*hgammas(iv) &
+ *(G31(ir,iv)*hpxis(k)*hgammas(m) &
+ +G33(ir,iv)*hxis(k)*hpgammas(m))
+
+ enddo
+ enddo
+
+ enddo
+ enddo
+
+ end subroutine compute_arrays_source
+
Modified: seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90 2007-12-07 23:46:47 UTC (rev 8444)
+++ seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90 2007-12-07 23:46:47 UTC (rev 8445)
@@ -1,19 +1,15 @@
-!=====================================================================
+
+!========================================================================
!
-! S p e c f e m 3 D B a s i n V e r s i o n 1 . 2
-! --------------------------------------------------
+! S P E C F E M 2 D Version 5.1
+! ------------------------------
!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology July 2004
+! Dimitri Komatitsch
+! Universite de Pau et des Pays de l'Adour, France
!
-! A signed non-commercial agreement is required to use this program.
-! Please check http://www.gps.caltech.edu/research/jtromp for details.
-! Free for non-commercial academic research ONLY.
-! This program is distributed WITHOUT ANY WARRANTY whatsoever.
-! Do not redistribute this program without written permission.
+! (c) January 2005
!
-!=====================================================================
+!========================================================================
program convolve_source_time_function
Modified: seismo/2D/SPECFEM2D/trunk/defarrays.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/defarrays.f90 2007-12-07 23:46:47 UTC (rev 8444)
+++ seismo/2D/SPECFEM2D/trunk/defarrays.f90 2007-12-07 23:46:47 UTC (rev 8445)
@@ -12,10 +12,9 @@
!========================================================================
subroutine defarrays(vpext,vsext,rhoext,density,elastcoef, &
- xigll,zigll,xix,xiz,gammax,gammaz,a11,a12, &
ibool,kmato,coord,npoin,rsizemin,rsizemax, &
cpoverdxmin,cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax, &
- vpmin,vpmax,readmodel,nspec,numat,source_type,ix_source,iz_source,ispec_source)
+ vpmin,vpmax,readmodel,nspec,numat)
! define all the arrays for the variational formulation
@@ -24,23 +23,13 @@
include "constants.h"
integer i,j,ispec,material,ipointnum,npoin,nspec,numat
- integer ix_source,iz_source,ispec_source,ir,is,source_type
integer kmato(nspec),ibool(NGLLX,NGLLX,nspec)
- double precision xix(NGLLX,NGLLZ,nspec)
- double precision xiz(NGLLX,NGLLZ,nspec)
- double precision gammax(NGLLX,NGLLZ,nspec)
- double precision gammaz(NGLLX,NGLLZ,nspec)
-
double precision density(numat),elastcoef(4,numat)
double precision coord(NDIM,npoin)
- double precision a11(NGLLX,NGLLX),a12(NGLLX,NGLLX)
-
- double precision xigll(NGLLX),zigll(NGLLZ)
-
double precision vpext(npoin)
double precision vsext(npoin)
double precision rhoext(npoin)
@@ -50,24 +39,18 @@
double precision kappa,cploc,csloc,x0,z0
double precision x1,z1,x2,z2,rdist1,rdist2,rapportmin,rapportmax
double precision lambdamin,lambdamax
- double precision flagxprime,flagzprime,sig0
double precision rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax,vpmin,vpmax
logical readmodel
- double precision, external :: lagrange_deriv_GLL
-
!
!-----------------------------------------------------------------------
!
!---- compute parameters for the spectral elements
- a11 = zero
- a12 = zero
-
vpmin = HUGEVAL
vsmin = HUGEVAL
vpmax = -HUGEVAL
@@ -181,29 +164,5 @@
print *,'********'
print *
-! seulement si source explosive
- if(source_type == 2) then
-
- if(ix_source == 1 .or. ix_source == NGLLX .or. iz_source == 1 .or. iz_source == NGLLX) &
- stop 'Explosive source on element edge'
-
-!---- definir a11 et a12 - dirac (schema en croix)
-
- sig0 = one
-
- do ir=1,NGLLX
- flagxprime = lagrange_deriv_GLL(ir-1,ix_source-1,xigll,NGLLX)
- a11(ir,iz_source) = a11(ir,iz_source) + sig0*xix(ix_source,iz_source,ispec_source)*flagxprime
- a12(ir,iz_source) = a12(ir,iz_source) + sig0*xiz(ix_source,iz_source,ispec_source)*flagxprime
- enddo
-
- do is=1,NGLLZ
- flagzprime = lagrange_deriv_GLL(is-1,iz_source-1,zigll,NGLLZ)
- a11(ix_source,is) = a11(ix_source,is) + sig0*gammax(ix_source,iz_source,ispec_source)*flagzprime
- a12(ix_source,is) = a12(ix_source,is) + sig0*gammaz(ix_source,iz_source,ispec_source)*flagzprime
- enddo
-
- endif
-
end subroutine defarrays
Modified: seismo/2D/SPECFEM2D/trunk/locate_receivers.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_receivers.f90 2007-12-07 23:46:47 UTC (rev 8444)
+++ seismo/2D/SPECFEM2D/trunk/locate_receivers.f90 2007-12-07 23:46:47 UTC (rev 8445)
@@ -171,7 +171,6 @@
write(IOUT,*) ' original x: ',sngl(st_xval(irec))
write(IOUT,*) ' original z: ',sngl(st_zval(irec))
write(IOUT,*) ' distance from source: ',sngl(distance_receiver)
-
write(IOUT,*) 'closest estimate found: ',sngl(final_distance(irec)),' m away'
write(IOUT,*) ' in element ',ispec_selected_rec(irec)
write(IOUT,*) ' at xi,gamma coordinates = ',xi_receiver(irec),gamma_receiver(irec)
Added: seismo/2D/SPECFEM2D/trunk/locate_source_force.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_source_force.f90 2007-12-07 23:46:47 UTC (rev 8444)
+++ seismo/2D/SPECFEM2D/trunk/locate_source_force.f90 2007-12-07 23:46:47 UTC (rev 8445)
@@ -0,0 +1,92 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.1
+! ------------------------------
+!
+! Dimitri Komatitsch
+! Universite de Pau et des Pays de l'Adour, France
+!
+! (c) January 2005
+!
+!========================================================================
+
+ subroutine locate_source_force(coord,ibool,npoin,nspec,x_source,z_source,source_type,ix_source,iz_source,ispec_source,iglob_source)
+
+!
+!----- calculer la position reelle de la source
+!
+
+ implicit none
+
+ include "constants.h"
+
+ integer npoin,nspec,source_type
+ integer ibool(NGLLX,NGLLZ,nspec)
+
+ double precision x_source,z_source
+ double precision coord(NDIM,npoin)
+
+ integer ip,ix,iz,numelem,ilowx,ilowz,ihighx,ihighz,ix_source,iz_source,ispec_source,iglob_source
+
+ double precision distminmax,distmin,xp,zp,dist
+
+ write(iout,200)
+
+ distminmax = -HUGEVAL
+
+ distmin = +HUGEVAL
+
+ ilowx = 1
+ ilowz = 1
+ 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
+
+! recherche du point de grille le plus proche
+ do numelem=1,nspec
+ do ix=ilowx,ihighx
+ do iz=ilowz,ihighz
+
+! numero global du point
+ ip=ibool(ix,iz,numelem)
+
+! coordonnees du point de grille
+ 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
+ if(dist < distmin) then
+ distmin = dist
+ iglob_source = ip
+ ix_source = ix
+ iz_source = iz
+ ispec_source = numelem
+ endif
+
+ enddo
+ enddo
+ enddo
+
+ distminmax = max(distmin,distminmax)
+
+ write(iout,150) x_source,z_source,coord(1,iglob_source),coord(2,iglob_source),distmin
+ write(iout,160) distminmax
+
+ 150 format(1x,f12.3,1x,f12.3,1x,f12.3,1x,f12.3,f12.3)
+ 160 format(/2x,'Maximum distance between asked and real =',f12.3)
+ 200 format(//1x,48('=')/,' = S o u r c e s ', &
+ 'r e a l p o s i t i o n s ='/1x,48('=')// &
+ ' Source x-asked z-asked x-obtain z-obtain dist'/)
+
+ end subroutine locate_source_force
+
Added: seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.f90 2007-12-07 23:46:47 UTC (rev 8444)
+++ seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.f90 2007-12-07 23:46:47 UTC (rev 8445)
@@ -0,0 +1,152 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.1
+! ------------------------------
+!
+! Dimitri Komatitsch
+! Universite de Pau et des Pays de l'Adour, France
+!
+! (c) January 2005
+!
+!========================================================================
+
+!----
+!---- locate_source_moment_tensor finds the correct position of the moment-tensor source
+!----
+
+ subroutine locate_source_moment_tensor(ibool,coord,nspec,npoin,xigll,zigll,x_source,z_source, &
+ ispec_selected_source,xi_source,gamma_source,coorg,knods,ngnod,npgeo)
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspec,npoin,ngnod,npgeo
+
+ integer knods(ngnod,nspec)
+ double precision coorg(NDIM,npgeo)
+
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
+
+! array containing coordinates of the points
+ double precision coord(NDIM,npoin)
+
+ integer i,j,ispec,iglob,iter_loop,ix_initial_guess,iz_initial_guess
+
+ double precision x_source,z_source,dist
+ double precision xi,gamma,dx,dz,dxi,dgamma
+
+! Gauss-Lobatto-Legendre points of integration
+ double precision xigll(NGLLX)
+ double precision zigll(NGLLZ)
+
+ double precision x,z,xix,xiz,gammax,gammaz,jacobian
+ double precision distmin,final_distance
+
+! source information
+ integer ispec_selected_source
+ double precision xi_source,gamma_source
+
+! **************
+
+ write(IOUT,*)
+ write(IOUT,*) '*******************************'
+ write(IOUT,*) ' locating moment-tensor source'
+ write(IOUT,*) '*******************************'
+ write(IOUT,*)
+
+! set distance to huge initial value
+ distmin=HUGEVAL
+
+ 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
+
+ iglob = ibool(i,j,ispec)
+ dist = sqrt((x_source-dble(coord(1,iglob)))**2 + (z_source-dble(coord(2,iglob)))**2)
+
+! keep this point if it is closer to the source
+ if(dist < distmin) then
+ distmin = dist
+ ispec_selected_source = ispec
+ ix_initial_guess = i
+ iz_initial_guess = j
+ endif
+
+ enddo
+ enddo
+
+! end of loop on all the spectral elements
+ enddo
+
+! ****************************************
+! find the best (xi,gamma) for each source
+! ****************************************
+
+! use initial guess in xi and gamma
+ xi = xigll(ix_initial_guess)
+ gamma = zigll(iz_initial_guess)
+
+! iterate to solve the non linear system
+ do iter_loop = 1,NUM_ITER
+
+! recompute jacobian for the new point
+ call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian,coorg,knods,ispec_selected_source,ngnod,nspec,npgeo)
+
+! compute distance to target location
+ dx = - (x - x_source)
+ dz = - (z - z_source)
+
+! compute increments
+ dxi = xix*dx + xiz*dz
+ dgamma = gammax*dx + gammaz*dz
+
+! update values
+ xi = xi + dxi
+ gamma = gamma + dgamma
+
+! impose that we stay in that element
+! (useful if user gives a source outside the mesh for instance)
+! we can go slightly outside the [1,1] segment since with finite elements
+! the polynomial solution is defined everywhere
+! this can be useful for convergence of itertive scheme with distorted elements
+ if (xi > 1.10d0) xi = 1.10d0
+ if (xi < -1.10d0) xi = -1.10d0
+ if (gamma > 1.10d0) gamma = 1.10d0
+ if (gamma < -1.10d0) gamma = -1.10d0
+
+! end of non linear iterations
+ enddo
+
+! compute final coordinates of point found
+ call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian,coorg,knods,ispec_selected_source,ngnod,nspec,npgeo)
+
+! store xi,gamma of point found
+ xi_source = xi
+ gamma_source = gamma
+
+! compute final distance between asked and found
+ final_distance = sqrt((x_source-x)**2 + (z_source-z)**2)
+
+ write(IOUT,*)
+ write(IOUT,*) 'Moment-tensor source:'
+
+ if(final_distance == HUGEVAL) stop 'error locating moment-tensor source'
+
+ write(IOUT,*) ' original x: ',sngl(x_source)
+ write(IOUT,*) ' original z: ',sngl(z_source)
+ write(IOUT,*) 'closest estimate found: ',sngl(final_distance),' m away'
+ write(IOUT,*) ' in element ',ispec_selected_source
+ write(IOUT,*) ' at xi,gamma coordinates = ',xi_source,gamma_source
+ write(IOUT,*)
+
+ write(IOUT,*)
+ write(IOUT,*) 'end of moment-tensor source detection'
+ write(IOUT,*)
+
+ end subroutine locate_source_moment_tensor
+
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.f90 2007-12-07 23:46:47 UTC (rev 8444)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.f90 2007-12-07 23:46:47 UTC (rev 8445)
@@ -44,7 +44,7 @@
! for the source
integer source_type,time_function_type
- double precision xs,zs,f0,t0,angle,factor
+ double precision xs,zs,f0,t0,angleforce,Mxx,Mzz,Mxz,factor
! arrays for the receivers
double precision, dimension(:), allocatable :: xrec,zrec
@@ -192,7 +192,10 @@
call read_value_integer(IIN_PAR,IGNORE_JUNK,source_type)
call read_value_integer(IIN_PAR,IGNORE_JUNK,time_function_type)
call read_value_double_precision(IIN_PAR,IGNORE_JUNK,f0)
- call read_value_double_precision(IIN_PAR,IGNORE_JUNK,angle)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,angleforce)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,Mxx)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,Mzz)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,Mxz)
call read_value_double_precision(IIN_PAR,IGNORE_JUNK,factor)
! if Dirac source time function, use a very thin Gaussian instead
@@ -207,7 +210,10 @@
print *,'Frequency, delay = ',f0,t0
print *,'Source type (1=force, 2=explosion): ',source_type
print *,'Time function type (1=Ricker, 2=First derivative, 3=Gaussian, 4=Dirac): ',time_function_type
- print *,'Angle of the source if force = ',angle
+ print *,'Angle of the source if force = ',angleforce
+ print *,'Mxx of the source if moment tensor = ',Mxx
+ print *,'Mzz of the source if moment tensor = ',Mzz
+ print *,'Mxz of the source if moment tensor = ',Mxz
print *,'Multiplying factor = ',factor
! read receivers line parameters
@@ -620,7 +626,7 @@
write(15,*) nt,dt
write(15,*) 'source'
- write(15,*) source_type,time_function_type,xs-xmin,zs,f0,t0,factor,angle
+ write(15,*) source_type,time_function_type,xs-xmin,zs,f0,t0,factor,angleforce,Mxx,Mzz,Mxz
write(15,*) 'Coordinates of macrobloc mesh (coorg):'
do j=0,nz
Deleted: seismo/2D/SPECFEM2D/trunk/positsource.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/positsource.f90 2007-12-07 23:46:47 UTC (rev 8444)
+++ seismo/2D/SPECFEM2D/trunk/positsource.f90 2007-12-07 23:46:47 UTC (rev 8445)
@@ -1,92 +0,0 @@
-
-!========================================================================
-!
-! S P E C F E M 2 D Version 5.1
-! ------------------------------
-!
-! Dimitri Komatitsch
-! Universite de Pau et des Pays de l'Adour, France
-!
-! (c) January 2005
-!
-!========================================================================
-
- subroutine positsource(coord,ibool,npoin,nspec,x_source,z_source,source_type,ix_source,iz_source,ispec_source,iglob_source)
-
-!
-!----- calculer la position reelle de la source
-!
-
- implicit none
-
- include "constants.h"
-
- integer npoin,nspec,source_type
- integer ibool(NGLLX,NGLLZ,nspec)
-
- double precision x_source,z_source
- double precision coord(NDIM,npoin)
-
- integer ip,ix,iz,numelem,ilowx,ilowz,ihighx,ihighz,ix_source,iz_source,ispec_source,iglob_source
-
- double precision distminmax,distmin,xp,zp,dist
-
- write(iout,200)
-
- distminmax = -HUGEVAL
-
- distmin = +HUGEVAL
-
- ilowx = 1
- ilowz = 1
- 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
-
-! recherche du point de grille le plus proche
- do numelem=1,nspec
- do ix=ilowx,ihighx
- do iz=ilowz,ihighz
-
-! numero global du point
- ip=ibool(ix,iz,numelem)
-
-! coordonnees du point de grille
- 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
- if(dist < distmin) then
- distmin = dist
- iglob_source = ip
- ix_source = ix
- iz_source = iz
- ispec_source = numelem
- endif
-
- enddo
- enddo
- enddo
-
- distminmax = max(distmin,distminmax)
-
- write(iout,150) x_source,z_source,coord(1,iglob_source),coord(2,iglob_source),distmin
- write(iout,160) distminmax
-
- 150 format(1x,f12.3,1x,f12.3,1x,f12.3,1x,f12.3,f12.3)
- 160 format(/2x,'Maximum distance between asked and real =',f12.3)
- 200 format(//1x,48('=')/,' = S o u r c e s ', &
- 'r e a l p o s i t i o n s ='/1x,48('=')// &
- ' Source x-asked z-asked x-obtain z-obtain dist'/)
-
- end subroutine positsource
-
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.f90 2007-12-07 23:46:47 UTC (rev 8444)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.f90 2007-12-07 23:46:47 UTC (rev 8445)
@@ -23,6 +23,7 @@
! - Dirac and Gaussian time sources and corresponding convolution routine
! - option for acoustic medium instead of elastic
! - receivers at any location, not only grid points
+! - moment-tensor source at any location, not only a grid point
! - color PNM snapshots
! - more flexible DATA/Par_file with any number of comment lines
! - Xsu scripts for seismograms
@@ -52,7 +53,8 @@
character(len=80) datlin
integer source_type,time_function_type
- double precision x_source,z_source,f0,t0,factor,angleforce
+ double precision x_source,z_source,xi_source,gamma_source,Mxx,Mzz,Mxz,f0,t0,factor,angleforce
+ double precision, dimension(NDIM,NGLLX,NGLLZ) :: sourcearray
double precision, dimension(:,:), allocatable :: coorg
double precision, dimension(:), allocatable :: coorgread
@@ -110,11 +112,8 @@
double precision, dimension(:), allocatable :: rmass,density,vpext,vsext,rhoext,displread,velocread,accelread
- double precision, dimension(:,:,:), allocatable :: shape2D,shape2D_display, &
- xix,xiz,gammax,gammaz,jacobian
+ double precision, dimension(:,:,:), allocatable :: shape2D,shape2D_display,xix,xiz,gammax,gammaz,jacobian
- double precision, dimension(:,:), allocatable :: a11,a12
-
double precision, dimension(:,:,:,:), allocatable :: dershape2D,dershape2D_display
integer, dimension(:,:,:), allocatable :: ibool
@@ -123,12 +122,12 @@
integer ie,k
- integer ispec_source,iglob_source,ix_source,iz_source
+ integer ispec_selected_source,iglob_source,ix_source,iz_source
double precision a,displnorm_all
double precision, dimension(:), allocatable :: source_time_function
double precision rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
- lambdal_Smin,lambdal_Smax,lambdal_Pmin,lambdal_Pmax,vpmin,vpmax
+ lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax,vpmin,vpmax
integer colors,numbers,subsamp,vecttype,itaff,nrec,sismostype
integer numat,ngnod,nspec,iptsdisp,nelemabs,nelemsurface
@@ -265,7 +264,7 @@
!---- read source information
!
read(IIN,40) datlin
- read(IIN,*) source_type,time_function_type,x_source,z_source,f0,t0,factor,angleforce
+ read(IIN,*) source_type,time_function_type,x_source,z_source,f0,t0,factor,angleforce,Mxx,Mzz,Mxz
!
!----- check the input
@@ -274,7 +273,7 @@
if (source_type == 1) then
write(IOUT,212) x_source,z_source,f0,t0,factor,angleforce
else if(source_type == 2) then
- write(IOUT,222) x_source,z_source,f0,t0,factor
+ write(IOUT,222) x_source,z_source,f0,t0,factor,Mxx,Mzz,Mxz
else
stop 'Unknown source type number !'
endif
@@ -319,8 +318,6 @@
allocate(gammax(NGLLX,NGLLZ,nspec))
allocate(gammaz(NGLLX,NGLLZ,nspec))
allocate(jacobian(NGLLX,NGLLZ,nspec))
- allocate(a11(NGLLX,NGLLZ))
- allocate(a12(NGLLX,NGLLZ))
allocate(flagrange(NGLLX,iptsdisp))
allocate(xinterp(iptsdisp,iptsdisp))
allocate(zinterp(iptsdisp,iptsdisp))
@@ -588,11 +585,25 @@
deltatover2 = HALF*deltat
deltatsquareover2 = HALF*deltat*deltat
-!
!---- definir la position reelle des points source et recepteurs
-!
- call positsource(coord,ibool,npoin,nspec,x_source,z_source,source_type,ix_source,iz_source,ispec_source,iglob_source)
+ if(source_type == 1) then
+! collocated force source
+ call locate_source_force(coord,ibool,npoin,nspec,x_source,z_source,source_type,ix_source,iz_source,ispec_selected_source,iglob_source)
+ else if(source_type == 2) then
+! moment-tensor source
+ call locate_source_moment_tensor(ibool,coord,nspec,npoin,xigll,zigll,x_source,z_source, &
+ ispec_selected_source,xi_source,gamma_source,coorg,knods,ngnod,npgeo)
+
+! compute source array for moment-tensor source
+ call compute_arrays_source(ispec_selected_source,xi_source,gamma_source,sourcearray, &
+ Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,nspec)
+
+ else
+ stop 'incorrect source type'
+ endif
+
+
! locate receivers in the mesh
call locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,st_xval,st_zval,ispec_selected_rec, &
xi_receiver,gamma_receiver,station_name,network_name,x_source,z_source,coorg,knods,ngnod,npgeo)
@@ -629,10 +640,9 @@
!---- define all arrays
!
call defarrays(vpext,vsext,rhoext,density,elastcoef, &
- xigll,zigll,xix,xiz,gammax,gammaz,a11,a12, &
ibool,kmato,coord,npoin,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,lambdal_Smin,lambdal_Smax,lambdal_Pmin,lambdal_Pmax, &
- vpmin,vpmax,readmodel,nspec,numat,source_type,ix_source,iz_source,ispec_source)
+ cpoverdxmin,cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax, &
+ vpmin,vpmax,readmodel,nspec,numat)
! build the global mass matrix once and for all
rmass(:) = ZERO
@@ -663,12 +673,10 @@
! convertir angle recepteurs en radians
anglerec = anglerec * pi / 180.d0
-!
!---- verifier le maillage, la stabilite et le nb de points par lambda
!---- seulement si la source en temps n'est pas un Dirac (sinon spectre non defini)
-!
if(time_function_type /= 4) call checkgrid(deltat,f0,t0,initialfield, &
- rsizemin,rsizemax,cpoverdxmin,cpoverdxmax,lambdal_Smin,lambdal_Smax,lambdal_Pmin,lambdal_Pmax)
+ rsizemin,rsizemax,cpoverdxmin,cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax)
!
!---- for color PNM images
@@ -1379,15 +1387,17 @@
accel(1,iglob_source) = accel(1,iglob_source) + source_time_function(it)
endif
-! explosion
+! moment tensor
else if(source_type == 2) then
- do i=1,NGLLX
- do j=1,NGLLX
- iglob = ibool(i,j,ispec_source)
- accel(1,iglob) = accel(1,iglob) + a11(i,j)*source_time_function(it)
- accel(2,iglob) = accel(2,iglob) + a12(i,j)*source_time_function(it)
+
+! add source array
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source)
+ accel(:,iglob) = accel(:,iglob) + sourcearray(:,i,j)*source_time_function(it)
enddo
enddo
+
endif
else
@@ -1743,10 +1753,10 @@
111 format('Sauvegarde deplacement temps t = ',1pe10.4,' s')
400 format(/1x,41('=')/,' = T i m e e v o l u t i o n l o o p ='/1x,41('=')/)
- 200 format(//1x,'C o n t r o l',/1x,34('='),//5x,&
+ 200 format(//1x,'C o n t r o l',/1x,13('='),//5x,&
'Number of spectral element control nodes. . .(npgeo) =',i8/5x, &
- 'Number of space dimensions . . . . . . . . . (NDIM) =',i8)
- 600 format(//1x,'C o n t r o l',/1x,34('='),//5x, &
+ 'Number of space dimensions. . . . . . . . . . (NDIM) =',i8)
+ 600 format(//1x,'C o n t r o l',/1x,13('='),//5x, &
'Display frequency . . . . . . . . . . . . . (itaff) = ',i5/ 5x, &
'Color display . . . . . . . . . . . . . . . (colors) = ',i5/ 5x, &
' == 0 black and white display ', / 5x, &
@@ -1754,24 +1764,24 @@
'Numbered mesh . . . . . . . . . . . . . . .(numbers) = ',i5/ 5x, &
' == 0 do not number the mesh ', /5x, &
' == 1 number the mesh ')
- 700 format(//1x,'C o n t r o l',/1x,34('='),//5x, &
- 'Seismograms recording type. . . . . . .(sismostype) = ',i6/5x, &
+ 700 format(//1x,'C o n t r o l',/1x,13('='),//5x, &
+ 'Seismograms recording type . . . . . . .(sismostype) = ',i6/5x, &
'Angle for first line of receivers. . . . .(anglerec) = ',f6.2)
- 750 format(//1x,'C o n t r o l',/1x,34('='),//5x, &
+ 750 format(//1x,'C o n t r o l',/1x,13('='),//5x, &
'Read external initial field or not . .(initialfield) = ',l6/5x, &
'Read external velocity model or not . . .(readmodel) = ',l6/5x, &
'Elastic simulation or acoustic. . . . . . .(ELASTIC) = ',l6/5x, &
'Turn anisotropy on or off. . . .(TURN_ANISOTROPY_ON) = ',l6/5x, &
'Turn attenuation on or off. . .(TURN_ATTENUATION_ON) = ',l6/5x, &
'Save grid in external file or not. . . .(outputgrid) = ',l6)
- 800 format(//1x,'C o n t r o l',/1x,34('='),//5x, &
- 'Vector display type . . . . . . . . . . .(vecttype) = ',i6/5x, &
+ 800 format(//1x,'C o n t r o l',/1x,13('='),//5x, &
+ 'Vector display type. . . . . . . . . . . .(vecttype) = ',i6/5x, &
'Percentage of cut for vector plots. . . . .(cutvect) = ',f6.2/5x, &
- 'Subsampling for velocity model display . .(subsamp) = ',i6)
+ 'Subsampling for velocity model display. . .(subsamp) = ',i6)
- 703 format(//' I t e r a t i o n s '/1x,29('='),//5x, &
+ 703 format(//' I t e r a t i o n s '/1x,19('='),//5x, &
'Number of time iterations . . . . .(NSTEP) =',i8,/5x, &
- 'Time step increment . . . . . . . .(deltat) =',1pe15.6,/5x, &
+ 'Time step increment. . . . . . . .(deltat) =',1pe15.6,/5x, &
'Total simulation duration . . . . . (ttot) =',1pe15.6)
107 format(/5x,'--> Isoparametric Spectral Elements <--',//)
@@ -1786,7 +1796,7 @@
'Number of absorbing elements . . . .(nelemabs) =',i7)
212 format(//,5x, &
- 'Source Type. . . . . . . . . . . . . . = Collocated Force',/5x, &
+ 'Source Type. . . . . . . . . . . . . . = Collocated Force',/5x, &
'X-position (meters). . . . . . . . . . =',1pe20.10,/5x, &
'Y-position (meters). . . . . . . . . . =',1pe20.10,/5x, &
'Fundamental frequency (Hz) . . . . . . =',1pe20.10,/5x, &
@@ -1794,12 +1804,15 @@
'Multiplying factor . . . . . . . . . . =',1pe20.10,/5x, &
'Angle from vertical direction (deg). . =',1pe20.10,/5x)
222 format(//,5x, &
- 'Source Type. . . . . . . . . . . . . . = Explosion',/5x, &
+ 'Source Type. . . . . . . . . . . . . . = Moment-tensor',/5x, &
'X-position (meters). . . . . . . . . . =',1pe20.10,/5x, &
'Y-position (meters). . . . . . . . . . =',1pe20.10,/5x, &
'Fundamental frequency (Hz) . . . . . . =',1pe20.10,/5x, &
'Time delay (s) . . . . . . . . . . . . =',1pe20.10,/5x, &
- 'Multiplying factor . . . . . . . . . . =',1pe20.10,/5x)
+ 'Multiplying factor . . . . . . . . . . =',1pe20.10,/5x, &
+ 'Mxx. . . . . . . . . . . . . . . . . . =',1pe20.10,/5x, &
+ 'Mzz. . . . . . . . . . . . . . . . . . =',1pe20.10,/5x, &
+ 'Mxz. . . . . . . . . . . . . . . . . . =',1pe20.10)
end program specfem2D
Modified: seismo/2D/SPECFEM2D/trunk/write_seismograms.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/write_seismograms.f90 2007-12-07 23:46:47 UTC (rev 8444)
+++ seismo/2D/SPECFEM2D/trunk/write_seismograms.f90 2007-12-07 23:46:47 UTC (rev 8445)
@@ -1,19 +1,15 @@
-!=====================================================================
+
+!========================================================================
!
-! S p e c f e m 3 D B a s i n V e r s i o n 1 . 2
-! --------------------------------------------------
+! S P E C F E M 2 D Version 5.1
+! ------------------------------
!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology July 2004
+! Dimitri Komatitsch
+! Universite de Pau et des Pays de l'Adour, France
!
-! A signed non-commercial agreement is required to use this program.
-! Please check http://www.gps.caltech.edu/research/jtromp for details.
-! Free for non-commercial academic research ONLY.
-! This program is distributed WITHOUT ANY WARRANTY whatsoever.
-! Do not redistribute this program without written permission.
+! (c) January 2005
!
-!=====================================================================
+!========================================================================
! write seismograms to text files
More information about the cig-commits
mailing list