[cig-commits] r20676 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Sep 3 10:01:13 PDT 2012
Author: dkomati1
Date: 2012-09-03 10:01:13 -0700 (Mon, 03 Sep 2012)
New Revision: 20676
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_Bielak_conditions.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_initialfield.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/write_seismograms.F90
Log:
fixed the case of an incident plane P wave
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_Bielak_conditions.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_Bielak_conditions.f90 2012-09-03 14:57:47 UTC (rev 20675)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_Bielak_conditions.f90 2012-09-03 17:01:13 UTC (rev 20676)
@@ -68,7 +68,6 @@
double precision c_inc, c_refl, time_offset, f0
double precision, dimension(NDIM) :: A_plane, B_plane, C_plane
-
! get the coordinates of the mesh point
x = coord(1,iglob) - x0_source
z = z0_source - coord(2,iglob)
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-09-03 14:57:47 UTC (rev 20675)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-09-03 17:01:13 UTC (rev 20676)
@@ -1025,10 +1025,6 @@
cpl = sqrt((kappal + 4._CUSTOM_REAL*mul_unrelaxed_elastic/3._CUSTOM_REAL)/rhol)
csl = sqrt(mul_unrelaxed_elastic/rhol)
-!!! DK DK
- c_inc = csl
-!!! DK DK
-
!--- left absorbing boundary
if(codeabs(IEDGE4,ispecabs)) then
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_initialfield.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_initialfield.F90 2012-09-03 14:57:47 UTC (rev 20675)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_initialfield.F90 2012-09-03 17:01:13 UTC (rev 20676)
@@ -131,20 +131,24 @@
endif
! allow negative anglesource(1): incidence from the right side of the domain
+ ! anglesource has been converted from degrees to radians before
anglesource_abs=abs(anglesource(1))
if (anglesource_abs > pi/2.d0 .and. source_type(1) /= 3) &
call exit_MPI("incorrect anglesource: must have 0 <= anglesource < 90")
- ! only implemented for homogeneous media therefore only 1 material supported
+ ! only implemented for homogeneous media therefore only one material supported
numat_local = numat
if (numat /= 1) then
-! if (myrank == 0) write(IOUT,*) 'not possible to have several materials with a plane wave, using the first material'
-! numat_local = 1
if (myrank == 0) then
- print *, 'This is not a homogenous model while plane wave initial condition'
- print *, ' is given by analytical formulae for a homogenous model.'
- print *, 'Use at your own risk!'
+ write(IOUT,*)
+ write(IOUT,*) 'It is not possible to have several materials with a plane wave, thus using the first material.'
+ write(IOUT,*) 'This is not a homogenous model, it contains ',numat,' materials'
+ write(IOUT,*) 'but the plane wave initial and boundary fields'
+ write(IOUT,*) 'are computed by analytical formulas for a homogenous model.'
+ write(IOUT,*) 'Thus use at your own risk!'
+ write(IOUT,*)
endif
+ numat_local = 1
endif
mu = poroelastcoef(2,1,numat_local)
@@ -233,9 +237,10 @@
C_plane(1)=0.d0; C_plane(2)=0.d0
endif
- ! correct A_plane and B_plane according to incident direction
+ ! correct A_plane, B_plane and C_plane according to incident direction
if (anglesource(1) < 0.) then
- A_plane(1)=-A_plane(1); B_plane(1)=-B_plane(1)
+ A_plane(1)=-A_plane(1)
+ B_plane(1)=-B_plane(1)
C_plane(1)=-C_plane(1)
endif
@@ -264,16 +269,11 @@
! check if zs = zmax (free surface)
if (myrank == 0 .and. abs(z_source(1)-zmax) > SMALLVALTOL) then
- print *, 'It is easier to set zs in SOURCE = zmax in interfacefile to keep track of the initial wavefront'
+ print *, 'It is sometimes easier to set zs in SOURCE = zmax in interfacefile to keep track of the initial wavefront'
endif
- ! initialize the time offset to put the plane wave not too close to the irregularity on the free surface
- ! add -t0 to match with the actual traveltime of plane waves
- if (abs(anglesource(1))<1.d0*pi/180.d0 .and. source_type(1)/=3) then
- time_offset = -1.d0*(zmax-zmin)/2.d0/c_inc - t0
- else
- time_offset = 0.d0 - t0
- endif
+ ! add -t0 to match the actual (analytical) traveltime of plane waves
+ time_offset = 0.d0 - t0
! to correctly center the initial plane wave in the mesh
x0_source = x_source(1)
@@ -282,7 +282,7 @@
if (myrank == 0) then
write(IOUT,*)
write(IOUT,*) 'You can modify the location of the initial plane wave by changing xs and zs in DATA/SOURCE.'
- write(IOUT,*) ' for instance: xs=',x_source(1),' zs=',z_source(1), ' (zs must be the height of the free surface)'
+ write(IOUT,*) ' for instance: xs=',x_source(1),' zs=',z_source(1), ' (zs can/should be the height of the free surface)'
write(IOUT,*)
endif
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/write_seismograms.F90 2012-09-03 14:57:47 UTC (rev 20675)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/write_seismograms.F90 2012-09-03 17:01:13 UTC (rev 20676)
@@ -412,66 +412,5 @@
deallocate(buffer_binary)
-!----
- if ( myrank == 0 ) then
-
-! ligne de recepteurs pour Xsu
- open(unit=11,file='OUTPUT_FILES/receiver_line_Xsu_XWindow',status='unknown')
-
-! subtract t0 from seismograms to get correct zero time
- write(11,110) FACTORXSU,NSTEP,deltat,-t0,nrec
-
- do irec=1,nrec
- ! this format statement might now work for larger meshes
- !write(11,"(f12.5)") st_xval(irec)
- write(11,*) st_xval(irec)
- if(irec < nrec) write(11,*) ','
- enddo
-
- if(seismotype == 1) then
- write(11,*) '@title="Ux at displacement@component"@<@Ux_file_single.bin'
- else if(seismotype == 2) then
- write(11,*) '@title="Ux at velocity@component"@<@Ux_file_single.bin'
- else
- write(11,*) '@title="Ux at acceleration@component"@<@Ux_file_single.bin'
- endif
-
- close(11)
-
-! script de visualisation
- open(unit=11,file='OUTPUT_FILES/show_receiver_line_Xsu',status='unknown')
- write(11,"('#!/bin/csh')")
- write(11,*)
- write(11,*) '/bin/rm -f tempfile receiver_line_Xsu_postscript'
- write(11,*) '# concatener toutes les lignes'
- write(11,*) 'tr -d ''\012'' <receiver_line_Xsu_XWindow >tempfile'
- write(11,*) '# remettre fin de ligne'
- write(11,*) 'echo " " >> tempfile'
- write(11,*) '# supprimer espaces, changer arobas, dupliquer'
- write(11,120)
- write(11,*) '/bin/rm -f tempfile'
- write(11,*) '# copier fichier pour sortie postscript'
- write(11,130)
- write(11,*) '/bin/rm -f tempfile'
- write(11,*) 'echo ''rm -f uxpoly.ps uzpoly.ps'' > tempfile'
- write(11,*) 'cat tempfile receiver_line_Xsu_postscript > tempfile2'
- write(11,*) '/bin/mv -f tempfile2 receiver_line_Xsu_postscript'
- write(11,*) '/bin/rm -f tempfile'
- write(11,*) '# executer commande xsu'
- write(11,*) 'sh receiver_line_Xsu_XWindow'
- write(11,*) '/bin/rm -f tempfile tempfile2'
- close(11)
-
-endif
-
-! formats
- 110 format('xwigb at xcur=',f8.2,'@n1=',i6,'@d1=',f15.8,'@f1=',f15.8,'@label1="Time@(s)"@label2="x@(m)"@n2=',i6,'@x2=')
-
- 120 format('sed -e ''1,$s/ //g'' -e ''1,$s/@/ /g'' -e ''1,1p'' -e ''$,$s/Ux/Uz/g'' <tempfile > receiver_line_Xsu_XWindow')
-
- 130 format('sed -e ''1,$s/xwigb/pswigp/g'' ', &
- '-e ''1,$s/Ux_file_single.bin/Ux_file_single.bin > uxpoly.ps/g'' ', &
- '-e ''1,$s/Uz_file_single.bin/Uz_file_single.bin > uzpoly.ps/g'' receiver_line_Xsu_XWindow > receiver_line_Xsu_postscript')
-
end subroutine write_seismograms
More information about the CIG-COMMITS
mailing list