[cig-commits] r18165 - in seismo/3D/SPECFEM3D/trunk: . doc/paper_draft/images src/generate_databases src/meshfem3D src/shared src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Sat Apr 2 03:11:08 PDT 2011
Author: danielpeter
Date: 2011-04-02 03:11:08 -0700 (Sat, 02 Apr 2011)
New Revision: 18165
Added:
seismo/3D/SPECFEM3D/trunk/doc/paper_draft/images/mesh_inner_outer.jpg
seismo/3D/SPECFEM3D/trunk/doc/paper_draft/images/mesh_inner_outer2.jpg
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_adj_source_frechet.f90
Removed:
seismo/3D/SPECFEM3D/trunk/doc/paper_draft/images/mesh_inner_outer.pdf
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/config.h.in
Modified:
seismo/3D/SPECFEM3D/trunk/config.h.in
seismo/3D/SPECFEM3D/trunk/src/generate_databases/Makefile.in
seismo/3D/SPECFEM3D/trunk/src/shared/compute_arrays_source.f90
seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in
seismo/3D/SPECFEM3D/trunk/src/shared/detect_surface.f90
seismo/3D/SPECFEM3D/trunk/src/shared/save_header_file.f90
seismo/3D/SPECFEM3D/trunk/src/shared/write_c_binary.c
seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_noDev.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
Log:
adds faster file i/o for kernel simulation using absorbing boundaries; puts routine compute_adj_source_frechet() into new file compute_adj_source_frechet.f90; fixes a bug with reconstructed b_R_xx etc. arrays when no attenuation is used for kernel simulations
Modified: seismo/3D/SPECFEM3D/trunk/config.h.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/config.h.in 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/config.h.in 2011-04-02 10:11:08 UTC (rev 18165)
@@ -49,3 +49,6 @@
/* Define to 1 if `lex' declares `yytext' as a `char *' by default, not a
`char[]'. */
#undef YYTEXT_POINTER
+
+/* Uncomment and define to select optimized file i/o for regional simulations */
+#define USE_MAP_FUNCTION
Added: seismo/3D/SPECFEM3D/trunk/doc/paper_draft/images/mesh_inner_outer.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D/trunk/doc/paper_draft/images/mesh_inner_outer.jpg
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Deleted: seismo/3D/SPECFEM3D/trunk/doc/paper_draft/images/mesh_inner_outer.pdf
===================================================================
(Binary files differ)
Added: seismo/3D/SPECFEM3D/trunk/doc/paper_draft/images/mesh_inner_outer2.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D/trunk/doc/paper_draft/images/mesh_inner_outer2.jpg
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/Makefile.in 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/Makefile.in 2011-04-02 10:11:08 UTC (rev 18165)
@@ -356,6 +356,6 @@
$O/param_reader.o: ${SHARED}/param_reader.c
${CC} -c $(CFLAGS) -o $O/param_reader.o ${SHARED}/param_reader.c -I../../
-$O/write_c_binary.o: ${SHARED}/write_c_binary.c
- ${CC} -c $(CFLAGS) -o $O/write_c_binary.o ${SHARED}/write_c_binary.c -I../../
+#$O/write_c_binary.o: ${SHARED}/write_c_binary.c
+# ${CC} -c $(CFLAGS) -o $O/write_c_binary.o ${SHARED}/write_c_binary.c -I../../
Deleted: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/config.h.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/config.h.in 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/config.h.in 2011-04-02 10:11:08 UTC (rev 18165)
@@ -1,34 +0,0 @@
-/* config.h.in. Generated from configure.ac by autoheader. */
-
-/* Define to dummy `main' function (if any) required to link to the Fortran
- libraries. */
-#undef FC_DUMMY_MAIN
-
-/* Define if F77 and FC dummy `main' functions are identical. */
-#undef FC_DUMMY_MAIN_EQ_F77
-
-/* Define to a macro mangling the given C identifier (in lower and upper
- case), which must not contain underscores, for linking with Fortran. */
-#undef FC_FUNC
-
-/* As FC_FUNC, but for C identifiers containing underscores. */
-#undef FC_FUNC_
-
-/* Define to alternate name for `main' routine that is called from a `main' in
- the Fortran libraries. */
-#undef FC_MAIN
-
-/* Define to the address where bug reports for this package should be sent. */
-#undef PACKAGE_BUGREPORT
-
-/* Define to the full name of this package. */
-#undef PACKAGE_NAME
-
-/* Define to the full name and version of this package. */
-#undef PACKAGE_STRING
-
-/* Define to the one symbol short name of this package. */
-#undef PACKAGE_TARNAME
-
-/* Define to the version of this package. */
-#undef PACKAGE_VERSION
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/compute_arrays_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/compute_arrays_source.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/compute_arrays_source.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -259,190 +259,8 @@
end subroutine compute_arrays_adjoint_source
-
! =======================================================================
-! compute the integrated derivatives of source parameters (M_jk and X_s)
-subroutine compute_adj_source_frechet(displ_s,Mxx,Myy,Mzz,Mxy,Mxz,Myz,eps_s,eps_m_s, &
- hxir,hetar,hgammar,hpxir,hpetar,hpgammar, hprime_xx,hprime_yy,hprime_zz, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
-
- implicit none
-
- include 'constants.h'
-
- ! input
- real(kind=CUSTOM_REAL) :: displ_s(NDIM,NGLLX,NGLLY,NGLLZ)
- double precision :: Mxx, Myy, Mzz, Mxy, Mxz, Myz
- ! output
- real(kind=CUSTOM_REAL) :: eps_s(NDIM,NDIM), eps_m_s(NDIM)
-
- ! auxilliary
- double precision :: hxir(NGLLX), hetar(NGLLY), hgammar(NGLLZ), &
- hpxir(NGLLX),hpetar(NGLLY),hpgammar(NGLLZ)
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
- real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
- real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
-
-! local variables
- real(kind=CUSTOM_REAL) :: tempx1l,tempx2l,tempx3l, tempy1l,tempy2l,tempy3l, &
- tempz1l,tempz2l,tempz3l, hp1, hp2, hp3, &
- xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl, &
- duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl, &
- xix_s,xiy_s,xiz_s,etax_s,etay_s,etaz_s,gammax_s,gammay_s,gammaz_s, &
- hlagrange_xi, hlagrange_eta, hlagrange_gamma, hlagrange
-
- real(kind=CUSTOM_REAL) :: eps(NDIM,NDIM), eps_array(NDIM,NDIM,NGLLX,NGLLY,NGLLZ), &
- eps_m_array(NGLLX,NGLLY,NGLLZ)
-
- integer i,j,k,l
-
-
-! first compute the strain at all the GLL points of the source element
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
-
- tempx1l = 0._CUSTOM_REAL
- tempx2l = 0._CUSTOM_REAL
- tempx3l = 0._CUSTOM_REAL
-
- tempy1l = 0._CUSTOM_REAL
- tempy2l = 0._CUSTOM_REAL
- tempy3l = 0._CUSTOM_REAL
-
- tempz1l = 0._CUSTOM_REAL
- tempz2l = 0._CUSTOM_REAL
- tempz3l = 0._CUSTOM_REAL
-
- do l=1,NGLLX
- hp1 = hprime_xx(i,l)
- tempx1l = tempx1l + displ_s(1,l,j,k)*hp1
- tempy1l = tempy1l + displ_s(2,l,j,k)*hp1
- tempz1l = tempz1l + displ_s(3,l,j,k)*hp1
-
- hp2 = hprime_yy(j,l)
- tempx2l = tempx2l + displ_s(1,i,l,k)*hp2
- tempy2l = tempy2l + displ_s(2,i,l,k)*hp2
- tempz2l = tempz2l + displ_s(3,i,l,k)*hp2
-
- hp3 = hprime_zz(k,l)
- tempx3l = tempx3l + displ_s(1,i,j,l)*hp3
- tempy3l = tempy3l + displ_s(2,i,j,l)*hp3
- tempz3l = tempz3l + displ_s(3,i,j,l)*hp3
- enddo
-
-! dudx
- xixl = xix(i,j,k)
- xiyl = xiy(i,j,k)
- xizl = xiz(i,j,k)
- etaxl = etax(i,j,k)
- etayl = etay(i,j,k)
- etazl = etaz(i,j,k)
- gammaxl = gammax(i,j,k)
- gammayl = gammay(i,j,k)
- gammazl = gammaz(i,j,k)
-
- duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
- duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
- duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
-
- duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
- duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
- duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
-
- duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
- duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
- duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
-
-! strain eps_jk
- eps(1,1) = duxdxl
- eps(1,2) = (duxdyl + duydxl) / 2
- eps(1,3) = (duxdzl + duzdxl) / 2
- eps(2,2) = duydyl
- eps(2,3) = (duydzl + duzdyl) / 2
- eps(3,3) = duzdzl
- eps(2,1) = eps(1,2)
- eps(3,1) = eps(1,3)
- eps(3,2) = eps(2,3)
-
- eps_array(:,:,i,j,k) = eps(:,:)
-
-! Mjk eps_jk
- eps_m_array(i,j,k) = Mxx * eps(1,1) + Myy * eps(2,2) + Mzz * eps(3,3) + &
- 2 * (Mxy * eps(1,2) + Mxz * eps(1,3) + Myz * eps(2,3))
-
- enddo
- enddo
- enddo
-
- ! interpolate the strain eps_s(:,:) from eps_array(:,:,i,j,k)
- eps_s = 0.
- xix_s = 0.; xiy_s = 0.; xiz_s = 0.
- etax_s = 0.; etay_s = 0.; etaz_s = 0.
- gammax_s = 0.; gammay_s = 0.; gammaz_s = 0.
-
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
-
- hlagrange = hxir(i)*hetar(j)*hgammar(k)
-
- eps_s(1,1) = eps_s(1,1) + eps_array(1,1,i,j,k)*hlagrange
- eps_s(1,2) = eps_s(1,2) + eps_array(1,2,i,j,k)*hlagrange
- eps_s(1,3) = eps_s(1,3) + eps_array(1,3,i,j,k)*hlagrange
- eps_s(2,2) = eps_s(2,2) + eps_array(2,2,i,j,k)*hlagrange
- eps_s(2,3) = eps_s(2,3) + eps_array(2,3,i,j,k)*hlagrange
- eps_s(3,3) = eps_s(3,3) + eps_array(3,3,i,j,k)*hlagrange
-
- xix_s = xix_s + xix(i,j,k)*hlagrange
- xiy_s = xiy_s + xiy(i,j,k)*hlagrange
- xiz_s = xiz_s + xiz(i,j,k)*hlagrange
- etax_s = etax_s + etax(i,j,k)*hlagrange
- etay_s = etay_s + etay(i,j,k)*hlagrange
- etaz_s = etaz_s + etaz(i,j,k)*hlagrange
- gammax_s = gammax_s + gammax(i,j,k)*hlagrange
- gammay_s = gammay_s + gammay(i,j,k)*hlagrange
- gammaz_s = gammaz_s + gammaz(i,j,k)*hlagrange
-
- enddo
- enddo
- enddo
-
-! for completion purpose, not used in specfem3D.f90
- eps_s(2,1) = eps_s(1,2)
- eps_s(3,1) = eps_s(1,3)
- eps_s(3,2) = eps_s(2,3)
-
-! compute the gradient of M_jk * eps_jk, and then interpolate it
-
- eps_m_s = 0.
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
-
- hlagrange_xi = hpxir(i)*hetar(j)*hgammar(k)
- hlagrange_eta = hxir(i)*hpetar(j)*hgammar(k)
- hlagrange_gamma = hxir(i)*hetar(j)*hpgammar(k)
-
- eps_m_s(1) = eps_m_s(1) + eps_m_array(i,j,k) * (hlagrange_xi * xix_s &
- + hlagrange_eta * etax_s + hlagrange_gamma * gammax_s)
- eps_m_s(2) = eps_m_s(2) + eps_m_array(i,j,k) * (hlagrange_xi * xiy_s &
- + hlagrange_eta * etay_s + hlagrange_gamma * gammay_s)
- eps_m_s(3) = eps_m_s(3) + eps_m_array(i,j,k) * (hlagrange_xi * xiz_s &
- + hlagrange_eta * etaz_s + hlagrange_gamma * gammaz_s)
-
- enddo
- enddo
- enddo
-
-end subroutine compute_adj_source_frechet
-
-! =======================================================================
-
! compute array for acoustic source
subroutine compute_arrays_source_acoustic(xi_source,eta_source,gamma_source,&
sourcearray,xigll,yigll,zigll,factor_source)
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in 2011-04-02 10:11:08 UTC (rev 18165)
@@ -84,8 +84,8 @@
! I/O unit for source and receiver vtk file
integer, parameter :: IOVTK = 98
! I/O unit for absorbing boundary snapshots
- integer, parameter :: IOABS = 31
- integer, parameter :: IOABS_AC = 32
+! integer, parameter :: IOABS = 31
+! integer, parameter :: IOABS_AC = 32
! I/O unit for plotting source time function
integer, Parameter :: IOSTF = 71
! I/O unit for interface file
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/detect_surface.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/detect_surface.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/detect_surface.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -61,7 +61,7 @@
!local parameters
integer, dimension(:), allocatable :: valence_external_mesh
- integer :: ispec,i,j,k,ii,jj,kk,iglob,ier
+ integer :: ispec,i,j,k,iglob,ier
! detecting surface points/elements (based on valence check on NGLL points) for external mesh
allocate(valence_external_mesh(nglob),stat=ier)
@@ -107,30 +107,10 @@
) then
iglob = ibool(i,j,k,ispec)
if (valence_external_mesh(iglob) == 1) then
- ispec_is_surface_external_mesh(ispec) = .true.
-
- ! sets flags for all gll points on this face
- if (k == 1 .or. k == NGLLZ) then
- do jj = 1, NGLLY
- do ii = 1, NGLLX
- iglob_is_surface_external_mesh(ibool(ii,jj,k,ispec)) = .true.
- enddo
- enddo
- endif
- if (j == 1 .or. j == NGLLY) then
- do kk = 1, NGLLZ
- do ii = 1, NGLLX
- iglob_is_surface_external_mesh(ibool(ii,j,kk,ispec)) = .true.
- enddo
- enddo
- endif
- if (i == 1 .or. i == NGLLX) then
- do kk = 1, NGLLZ
- do jj = 1, NGLLY
- iglob_is_surface_external_mesh(ibool(i,jj,kk,ispec)) = .true.
- enddo
- enddo
- endif
+ ! sets surface flags for element and global points
+ call ds_set_surface_flags(nspec,ispec_is_surface_external_mesh, &
+ nglob,iglob_is_surface_external_mesh, &
+ i,j,k,ispec,ibool)
endif
endif
@@ -175,6 +155,60 @@
!-------------------------------------------------------------------------------------------------
!
+ subroutine ds_set_surface_flags(nspec,ispec_is_surface_external_mesh, &
+ nglob,iglob_is_surface_external_mesh, &
+ i,j,k,ispec,ibool)
+
+ ! put this into separate subroutine to compile faster, otherwise compilers will try to unroll all do loops
+
+ implicit none
+
+ include "constants.h"
+
+ ! global indexing
+ integer :: nglob,nspec
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec):: ibool
+ integer :: i,j,k,ispec
+
+ ! surface flags
+ logical, dimension(nspec) :: ispec_is_surface_external_mesh
+ logical, dimension(nglob) :: iglob_is_surface_external_mesh
+
+ ! local parameters
+ integer :: kk,jj,ii
+
+ ! sets surface flag for element
+ ispec_is_surface_external_mesh(ispec) = .true.
+
+ ! sets flags for all gll points on this face
+ if (k == 1 .or. k == NGLLZ) then
+ do jj = 1, NGLLY
+ do ii = 1, NGLLX
+ iglob_is_surface_external_mesh(ibool(ii,jj,k,ispec)) = .true.
+ enddo
+ enddo
+ endif
+ if (j == 1 .or. j == NGLLY) then
+ do kk = 1, NGLLZ
+ do ii = 1, NGLLX
+ iglob_is_surface_external_mesh(ibool(ii,j,kk,ispec)) = .true.
+ enddo
+ enddo
+ endif
+ if (i == 1 .or. i == NGLLX) then
+ do kk = 1, NGLLZ
+ do jj = 1, NGLLY
+ iglob_is_surface_external_mesh(ibool(i,jj,kk,ispec)) = .true.
+ enddo
+ enddo
+ endif
+
+ end subroutine ds_set_surface_flags
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine detect_surface_cross_section(NPROC,nglob,nspec,ibool,&
ispec_is_surface_external_mesh, &
iglob_is_surface_external_mesh, &
@@ -229,11 +263,11 @@
real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord_face,ycoord_face,zcoord_face
real(kind=CUSTOM_REAL) :: mindist,normal(NDIM)
integer, dimension(:), allocatable :: valence_external_mesh
- integer,dimension(3,NGLLX,NGLLX) :: face_ijk
- integer :: ispec,i,j,k,ii,jj,kk,iglob,ier,count
+
+ integer :: ispec,i,j,k,iglob,ier,count
integer :: iface,icorner
logical, dimension(:),allocatable :: ispec_has_points
- logical :: has_face
+
! corners indices of reference cube faces
integer,dimension(3,4),parameter :: iface1_corner_ijk = &
reshape((/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ /),(/3,4/)) ! xmin
@@ -332,56 +366,15 @@
! considers only points in same process or, if point is shared between two processes,
! only with higher process ranks than itself
- if (valence_external_mesh(iglob) == myrank+1 .or. valence_external_mesh(iglob) > 2*(myrank+1) ) then
-
- has_face = .false.
-
-
- ! sets flags for all gll points on a face and makes sure it's not inside the element
- ! zmin & zmax face
- if ((k == 1 .or. k == NGLLZ) .and. valence_external_mesh(ibool(3,3,k,ispec)) >= 1 ) then
- has_face = .true.
- do jj = 1, NGLLY
- do ii = 1, NGLLX
- iglob_is_surface_external_mesh(ibool(ii,jj,k,ispec)) = .true.
- ! resets valence to count face only once
- valence_external_mesh(ibool(ii,jj,k,ispec)) = -1
- enddo
- enddo
- endif
-
- ! ymin & ymax
- if ((j == 1 .or. j == NGLLY) .and. valence_external_mesh(ibool(3,j,3,ispec)) >= 1) then
- has_face = .true.
- do kk = 1, NGLLZ
- do ii = 1, NGLLX
- iglob_is_surface_external_mesh(ibool(ii,j,kk,ispec)) = .true.
- ! resets valence to count face only once
- valence_external_mesh(ibool(ii,j,kk,ispec)) = -1
- enddo
- enddo
- endif
-
- ! xmin & xmax
- if ((i == 1 .or. i == NGLLX) .and. valence_external_mesh(ibool(i,3,3,ispec)) >= 1) then
- has_face = .true.
- do kk = 1, NGLLZ
- do jj = 1, NGLLY
- iglob_is_surface_external_mesh(ibool(i,jj,kk,ispec)) = .true.
- ! resets valence to count face only once
- valence_external_mesh(ibool(i,jj,kk,ispec)) = -1
- enddo
- enddo
- endif
-
-
- ! sets flag for element
- if( has_face ) then
- ispec_is_surface_external_mesh(ispec) = .true.
- count = count+1
- endif
-
+ if (valence_external_mesh(iglob) == myrank+1 &
+ .or. valence_external_mesh(iglob) > 2*(myrank+1) ) then
+ ! sets surface flags for cross section
+ call ds_set_cross_section_flags(nspec,ispec_is_surface_external_mesh, &
+ nglob,iglob_is_surface_external_mesh, &
+ i,j,k,ispec,ibool, &
+ valence_external_mesh,count)
endif
+
endif
enddo
enddo
@@ -453,20 +446,11 @@
valence_external_mesh(ibool(i,j,k,ispec)) /= -1 ) then
! checks face normal points in similar direction as cross-section normal
if( abs(normal(1)) > 0.6 ) then
- call get_element_face_gll_indices(iface,face_ijk,NGLLX,NGLLX)
- do jj = 1, NGLLY
- do ii = 1, NGLLX
- i = face_ijk(1,ii,jj)
- j = face_ijk(2,ii,jj)
- k = face_ijk(3,ii,jj)
- ! sets iglob flag on face points
- iglob_is_surface_external_mesh(ibool(i,j,k,ispec)) = .true.
- ! sets ispec flag
- ispec_is_surface_external_mesh(ispec) = .true.
- ! resets valence
- valence_external_mesh(ibool(i,j,k,ispec)) = -1
- enddo
- enddo
+ ! sets surfaces flags
+ call ds_set_plane_flags(iface,ispec, &
+ nspec,ispec_is_surface_external_mesh, &
+ nglob,iglob_is_surface_external_mesh, &
+ ibool,valence_external_mesh)
endif
endif
@@ -480,20 +464,11 @@
valence_external_mesh(ibool(i,j,k,ispec)) /= -1) then
! checks face normal points in similar direction as cross-section normal
if( abs(normal(2)) > 0.6 ) then
- call get_element_face_gll_indices(iface,face_ijk,NGLLX,NGLLX)
- do jj = 1, NGLLY
- do ii = 1, NGLLX
- i = face_ijk(1,ii,jj)
- j = face_ijk(2,ii,jj)
- k = face_ijk(3,ii,jj)
- ! sets iglob flag on face points
- iglob_is_surface_external_mesh(ibool(i,j,k,ispec)) = .true.
- ! sets ispec flag
- ispec_is_surface_external_mesh(ispec) = .true.
- ! resets valence
- valence_external_mesh(ibool(i,j,k,ispec)) = -1
- enddo
- enddo
+ ! sets surfaces flags
+ call ds_set_plane_flags(iface,ispec, &
+ nspec,ispec_is_surface_external_mesh, &
+ nglob,iglob_is_surface_external_mesh, &
+ ibool,valence_external_mesh)
endif
endif
@@ -507,20 +482,11 @@
valence_external_mesh(ibool(i,j,k,ispec)) /= -1) then
! checks face normal points in similar direction as cross-section normal
if( abs(normal(3)) > 0.6 ) then
- call get_element_face_gll_indices(iface,face_ijk,NGLLX,NGLLX)
- do jj = 1, NGLLY
- do ii = 1, NGLLX
- i = face_ijk(1,ii,jj)
- j = face_ijk(2,ii,jj)
- k = face_ijk(3,ii,jj)
- ! sets iglob flag on face points
- iglob_is_surface_external_mesh(ibool(i,j,k,ispec)) = .true.
- ! sets ispec flag
- ispec_is_surface_external_mesh(ispec) = .true.
- ! resets valence
- valence_external_mesh(ibool(i,j,k,ispec)) = -1
- enddo
- enddo
+ ! sets surfaces flags
+ call ds_set_plane_flags(iface,ispec, &
+ nspec,ispec_is_surface_external_mesh, &
+ nglob,iglob_is_surface_external_mesh, &
+ ibool,valence_external_mesh)
endif
endif
@@ -562,10 +528,138 @@
end subroutine detect_surface_cross_section
+
!
!-------------------------------------------------------------------------------------------------
!
+ subroutine ds_set_cross_section_flags(nspec,ispec_is_surface_external_mesh, &
+ nglob,iglob_is_surface_external_mesh, &
+ i,j,k,ispec,ibool, &
+ valence_external_mesh,count)
+
+ ! put this into separate subroutine to compile faster, otherwise compilers will try to unroll all do loops
+
+ implicit none
+
+ include "constants.h"
+
+ ! global indexing
+ integer :: nglob,nspec
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec):: ibool
+ integer :: i,j,k,ispec,count
+
+ integer, dimension(nglob) :: valence_external_mesh
+
+ ! surface flags
+ logical, dimension(nspec) :: ispec_is_surface_external_mesh
+ logical, dimension(nglob) :: iglob_is_surface_external_mesh
+
+ ! local parameters
+ integer :: kk,jj,ii
+ logical :: has_face
+
+ ! initialize element flag
+ has_face = .false.
+
+ ! sets flags for all gll points on a face and makes sure it's not inside the element
+ ! zmin & zmax face
+ if ((k == 1 .or. k == NGLLZ) .and. valence_external_mesh(ibool(3,3,k,ispec)) >= 1 ) then
+ has_face = .true.
+ do jj = 1, NGLLY
+ do ii = 1, NGLLX
+ iglob_is_surface_external_mesh(ibool(ii,jj,k,ispec)) = .true.
+ ! resets valence to count face only once
+ valence_external_mesh(ibool(ii,jj,k,ispec)) = -1
+ enddo
+ enddo
+ endif
+
+ ! ymin & ymax
+ if ((j == 1 .or. j == NGLLY) .and. valence_external_mesh(ibool(3,j,3,ispec)) >= 1) then
+ has_face = .true.
+ do kk = 1, NGLLZ
+ do ii = 1, NGLLX
+ iglob_is_surface_external_mesh(ibool(ii,j,kk,ispec)) = .true.
+ ! resets valence to count face only once
+ valence_external_mesh(ibool(ii,j,kk,ispec)) = -1
+ enddo
+ enddo
+ endif
+
+ ! xmin & xmax
+ if ((i == 1 .or. i == NGLLX) .and. valence_external_mesh(ibool(i,3,3,ispec)) >= 1) then
+ has_face = .true.
+ do kk = 1, NGLLZ
+ do jj = 1, NGLLY
+ iglob_is_surface_external_mesh(ibool(i,jj,kk,ispec)) = .true.
+ ! resets valence to count face only once
+ valence_external_mesh(ibool(i,jj,kk,ispec)) = -1
+ enddo
+ enddo
+ endif
+
+ ! sets flag for element to indicate that it has a face on surface
+ if( has_face ) then
+ ispec_is_surface_external_mesh(ispec) = .true.
+ count = count+1
+ endif
+
+ end subroutine ds_set_cross_section_flags
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine ds_set_plane_flags(iface,ispec, &
+ nspec,ispec_is_surface_external_mesh, &
+ nglob,iglob_is_surface_external_mesh, &
+ ibool,valence_external_mesh)
+
+ ! put this into separate subroutine to compile faster, otherwise compilers will try to unroll all do loops
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: iface,ispec
+
+ ! global indexing
+ integer :: nglob,nspec
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec):: ibool
+ integer, dimension(nglob) :: valence_external_mesh
+
+ ! surface flags
+ logical, dimension(nspec) :: ispec_is_surface_external_mesh
+ logical, dimension(nglob) :: iglob_is_surface_external_mesh
+
+
+ ! local parameters
+ integer :: jj,ii,i,j,k
+ integer,dimension(3,NGLLX,NGLLX) :: face_ijk
+
+ call get_element_face_gll_indices(iface,face_ijk,NGLLX,NGLLX)
+
+ do jj = 1, NGLLY
+ do ii = 1, NGLLX
+ i = face_ijk(1,ii,jj)
+ j = face_ijk(2,ii,jj)
+ k = face_ijk(3,ii,jj)
+ ! sets iglob flag on face points
+ iglob_is_surface_external_mesh(ibool(i,j,k,ispec)) = .true.
+ ! sets ispec flag
+ ispec_is_surface_external_mesh(ispec) = .true.
+ ! resets valence
+ valence_external_mesh(ibool(i,j,k,ispec)) = -1
+ enddo
+ enddo
+
+ end subroutine ds_set_plane_flags
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine detect_surface_PNM_GIF_image(NPROC,nglob,nspec,ibool,&
ispec_is_image_surface, &
iglob_is_image_surface, &
@@ -676,7 +770,3 @@
num_iglob_image_surface = count
end subroutine detect_surface_PNM_GIF_image
-
-
-
-
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/save_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/save_header_file.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/save_header_file.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -106,11 +106,6 @@
write(IOUT,*) '! NSPEC_ADJOINT = ', 1
write(IOUT,*) '! NGLOB_ADJOINT = ', 1
endif
- if (ATTENUATION .and. SIMULATION_TYPE == 3) then
- write(IOUT,*) '! NSPEC_ATT_AND_KERNEL = ', NSPEC_AB
- else
- write(IOUT,*) '! NSPEC_ATT_AND_KERNEL = ', 1
- endif
write(IOUT,*) '! '
write(IOUT,*) '! approximate least memory needed by the solver:'
write(IOUT,*) '! ----------------------------------------------'
@@ -118,7 +113,7 @@
write(IOUT,*) '! size of static arrays for the biggest slice = ',static_memory_size/1048576.d0,' MB'
write(IOUT,*) '! = ',static_memory_size/1073741824.d0,' GB'
write(IOUT,*) '!'
- write(IOUT,*) '! (should be below and typically equal to 80% of 1.5 GB = 1.2 GB on pangu'
+ write(IOUT,*) '! (should be below to 80% of 1.5 GB = 1.2 GB on pangu'
write(IOUT,*) '! at Caltech, and below and typically equal to 85% of 2 GB = 1.7 GB'
write(IOUT,*) '! on Marenostrum in Barcelona)'
write(IOUT,*) '! (if significantly more, the job will not run by lack of memory)'
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/write_c_binary.c
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/write_c_binary.c 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/write_c_binary.c 2011-04-02 10:11:08 UTC (rev 18165)
@@ -64,3 +64,528 @@
int dummy_unused_variable = write(fd, z, sizeof(float));
}
+
+/* ---------------------------------------
+
+ IO performance test
+
+ Software Optimization for High Performance Computing: Creating Faster Applications
+
+ By Isom L. Crawford and Kevin R. Wadleigh
+ Jul 18, 2003
+
+ - uses functions fopen/fread/fwrite for binary file I/O
+
+ --------------------------------------- */
+
+#define __USE_GNU
+#include <string.h>
+#include <regex.h>
+
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+
+/* fastest performance on nehalem nodes:
+
+ Linux 2.6.18-164.11.1.el5 #1 SMP Wed Jan 20 10:04:55 EST 2010 x86_64 x86_64 x86_64 GNU/Linux
+
+ achieved with 16 KB buffers: */
+
+//#define MAX_B 65536 // 64 KB
+//#define MAX_B 32768 // 32 KB
+#define MAX_B 16384 // 16 KB
+//#define MAX_B 8192 // 8 KB
+
+// absorbing files: instead of passing file descriptor, we use the array index
+// index 0 for elastic domain file
+// index 1 for acoutic domain file
+// (reserved, but unused yet) index 2 - for NOISE_TOMOGRAPHY (SURFACE_MOVIE)
+#define ABS_FILEID 3 // number of absorbing files, note: in C we start indexing arrays from 0
+
+// file points
+static FILE * fp_abs[ABS_FILEID];
+// file work buffers
+static char * work_buffer[ABS_FILEID];
+
+
+//void
+//FC_FUNC_(open_file_abs_r_fbin,OPEN_FILE_ABS_R_FBIN)(int *fid, char *filename,int *length, int *filesize){
+void open_file_abs_r_fbin(int *fid, char *filename,int *length, int *filesize){
+
+ // opens file for read access
+
+ //This sequence assigns the MAX_B array work_buffer to the file pointer
+ // to be used for its buffering. performance should benefit.
+ char * fncopy;
+ char * blank;
+ FILE *ft;
+
+ // checks filesize
+ if( *filesize == 0 ){
+ perror("Error file size for reading");
+ exit(EXIT_FAILURE);
+ }
+
+ // Trim the file name.
+ fncopy = strndup(filename, *length);
+ blank = strchr(fncopy, ' ');
+ if (blank != NULL) {
+ fncopy[blank - fncopy] = '\0';
+ }
+
+ // opens file
+ ft = fopen( fncopy, "r+" );
+ if( ft == NULL ) { perror("fopen"); exit(-1); }
+
+ // sets mode for full buffering
+ work_buffer[*fid] = (char *)malloc(MAX_B);
+ setvbuf( ft, work_buffer[*fid], _IOFBF, (size_t)MAX_B );
+
+ // stores file index id fid: from 0 to 8
+ fp_abs[*fid] = ft;
+
+ free(fncopy);
+}
+
+//void
+//FC_FUNC_(open_file_abs_w_fbin,OPEN_FILE_ABS_W_FBIN)(int *fid, char *filename, int *length, int *filesize){
+void open_file_abs_w_fbin(int *fid, char *filename, int *length, int *filesize){
+
+ // opens file for write access
+
+ //This sequence assigns the MAX_B array work_buffer to the file pointer
+ // to be used for its buffering. performance should benefit.
+ char * fncopy;
+ char * blank;
+ FILE *ft;
+
+ // checks filesize
+ if( *filesize == 0 ){
+ perror("Error file size for reading");
+ exit(EXIT_FAILURE);
+ }
+
+ // Trim the file name.
+ fncopy = strndup(filename, *length);
+ blank = strchr(fncopy, ' ');
+ if (blank != NULL) {
+ fncopy[blank - fncopy] = '\0';
+ }
+
+ // opens file
+ ft = fopen( fncopy, "w+" );
+ if( ft == NULL ) { perror("fopen"); exit(-1); }
+
+ // sets mode for full buffering
+ work_buffer[*fid] = (char *)malloc(MAX_B);
+ setvbuf( ft, work_buffer[*fid], _IOFBF, (size_t)MAX_B );
+
+ // stores file index id fid: from 0 to 8
+ fp_abs[*fid] = ft;
+
+ free(fncopy);
+
+}
+
+//void
+//FC_FUNC_(close_file_abs_fbin,CLOSE_FILE_ABS_FBIN)(int * fid){
+void close_file_abs_fbin(int * fid){
+
+ // closes file
+
+ fclose(fp_abs[*fid]);
+
+ free(work_buffer[*fid]);
+
+}
+
+//void
+//FC_FUNC_(write_abs_fbin,WRITE_ABS_FBIN)(int *fid, void *buffer, int *length, int *index){
+void write_abs_fbin(int *fid, void *buffer, int *length, int *index){
+
+ // writes binary file data in chunks of MAX_B
+
+ FILE *ft;
+ int itemlen,remlen,donelen,ret;
+ void *buf;
+
+ // file pointer
+ ft = fp_abs[*fid];
+
+ donelen = 0;
+ remlen = *length;
+ buf = buffer;
+ ret = 0;
+
+ //float dat[2];
+ //memcpy(dat,buffer,*length);
+ //printf("buffer: %f %f\n",dat[0],dat[1]);
+
+ // writes items of maximum MAX_B to the file
+ while (remlen > 0){
+
+ itemlen = MIN(remlen,MAX_B);
+ ret = fwrite(buf,1,itemlen,ft);
+ if (ret > 0){
+ donelen = donelen + ret;
+ remlen = remlen - MAX_B;
+ buf += MAX_B;
+ }
+ else{
+ remlen = 0;
+ }
+ }
+
+}
+
+//void
+//FC_FUNC_(read_abs_fbin,READ_ABS_FBIN)(int *fid, void *buffer, int *length, int *index){
+void read_abs_fbin(int *fid, void *buffer, int *length, int *index){
+
+ // reads binary file data in chunks of MAX_B
+
+ FILE *ft;
+ int ret,itemlen,remlen,donelen,pos;
+ void *buf;
+
+ // file pointer
+ ft = fp_abs[*fid];
+
+ // positions file pointer (for reverse time access)
+ pos = (*length) * (*index -1 );
+ fseek(ft, pos , SEEK_SET);
+
+ donelen = 0;
+ remlen = *length;
+ buf = buffer;
+ ret = 0;
+
+ // reads items of maximum MAX_B to the file
+ while (remlen > 0){
+
+ // checks end of file
+ if (ferror(ft) || feof(ft)) return;
+
+ itemlen = MIN(remlen,MAX_B);
+ ret = fread(buf,1,itemlen,ft);
+
+ if (ferror(ft) || feof(ft)) return;
+
+ if (ret > 0){
+ donelen = donelen + ret;
+ remlen = remlen - MAX_B;
+ buf += MAX_B;
+ }
+ else{
+ remlen = 0;
+ }
+ }
+
+ //float dat[2];
+ //memcpy(dat,buffer,*length);
+ //printf("return buffer: %f %f\n",dat[0],dat[1]);
+}
+
+
+
+
+/* ---------------------------------------
+
+ IO performance test
+
+
+ A Performance Comparison of "read" and "mmap" in the Solaris 8 OS
+
+ By Oyetunde Fadele, September 2002
+
+ http://developers.sun.com/solaris/articles/read_mmap.html
+
+ or
+
+ High-performance network programming, Part 2: Speed up processing at both the client and server
+
+ by Girish Venkatachalam
+
+ http://www.ibm.com/developerworks/aix/library/au-highperform2/
+
+
+ - uses functions mmap/memcpy for mapping file I/O
+
+ ------------------------------------- */
+
+
+#include <errno.h>
+#include <limits.h>
+#include <sys/mman.h>
+
+// file maps
+static char * map_abs[ABS_FILEID];
+// file descriptors
+static int map_fd_abs[ABS_FILEID];
+// file sizes
+static int filesize_abs[ABS_FILEID];
+
+//void
+//FC_FUNC_(open_file_abs_w_map,OPEN_FILE_ABS_W_MAP)(int *fid, char *filename, int *length, int *filesize){
+void open_file_abs_w_map(int *fid, char *filename, int *length, int *filesize){
+
+ // opens file for write access
+
+ int ft;
+ int result;
+ char *map;
+ char *fncopy;
+ char *blank;
+
+ // checks filesize
+ if( *filesize == 0 ){
+ perror("Error file size for writing");
+ exit(EXIT_FAILURE);
+ }
+
+ // Trim the file name.
+ fncopy = strndup(filename, *length);
+ blank = strchr(fncopy, ' ');
+ if (blank != NULL) {
+ fncopy[blank - fncopy] = '\0';
+ }
+
+ /* Open a file for writing.
+ * - Creating the file if it doesn't exist.
+ * - Truncating it to 0 size if it already exists. (not really needed)
+ *
+ * Note: "O_WRONLY" mode is not sufficient when mmaping.
+ */
+ ft = open(fncopy, O_RDWR | O_CREAT | O_TRUNC, (mode_t)0600);
+ if (ft == -1) {
+ perror("Error opening file for writing");
+ exit(EXIT_FAILURE);
+ }
+
+ // file index id fid: from 0 to 8
+ map_fd_abs[*fid] = ft;
+
+ free(fncopy);
+
+
+ /* Stretch the file size to the size of the (mmapped) array of ints
+ */
+ filesize_abs[*fid] = *filesize;
+ result = lseek(ft, filesize_abs[*fid] - 1, SEEK_SET);
+ if (result == -1) {
+ close(ft);
+ perror("Error calling fseek() to 'stretch' the file");
+ exit(EXIT_FAILURE);
+ }
+
+ //printf("file length: %d \n",filesize_abs[*fid]);
+
+
+ /* Something needs to be written at the end of the file to
+ * have the file actually have the new size.
+ * Just writing an empty string at the current file position will do.
+ *
+ * Note:
+ * - The current position in the file is at the end of the stretched
+ * file due to the call to lseek().
+ * - An empty string is actually a single '\0' character, so a zero-byte
+ * will be written at the last byte of the file.
+ */
+ result = write(ft, "", 1);
+ if (result != 1) {
+ close(ft);
+ perror("Error writing last byte of the file");
+ exit(EXIT_FAILURE);
+ }
+
+ /* Now the file is ready to be mmapped.
+ */
+ map = mmap(0, filesize_abs[*fid], PROT_READ | PROT_WRITE, MAP_SHARED, ft, 0);
+ if (map == MAP_FAILED) {
+ close(ft);
+ perror("Error mmapping the file");
+ exit(EXIT_FAILURE);
+ }
+
+ map_abs[*fid] = map;
+
+ //printf("file map: %d\n",*fid);
+
+}
+
+//void
+//FC_FUNC_(open_file_abs_r_map,OPEN_FILE_ABS_R_MAP)(int *fid, char *filename,int *length, int *filesize){
+void open_file_abs_r_map(int *fid, char *filename,int *length, int *filesize){
+
+ // opens file for read access
+ char * fncopy;
+ char * blank;
+ int ft;
+ char *map;
+
+ // checks filesize
+ if( *filesize == 0 ){
+ perror("Error file size for reading");
+ exit(EXIT_FAILURE);
+ }
+
+ // Trim the file name.
+ fncopy = strndup(filename, *length);
+ blank = strchr(fncopy, ' ');
+ if (blank != NULL) {
+ fncopy[blank - fncopy] = '\0';
+ }
+
+
+ ft = open(fncopy, O_RDONLY);
+ if (ft == -1) {
+ perror("Error opening file for reading");
+ exit(EXIT_FAILURE);
+ }
+
+ // file index id fid: from 0 to 8
+ map_fd_abs[*fid] = ft;
+
+ free(fncopy);
+
+ filesize_abs[*fid] = *filesize;
+
+ map = mmap(0, filesize_abs[*fid], PROT_READ, MAP_SHARED, ft, 0);
+ if (map == MAP_FAILED) {
+ close(ft);
+ perror("Error mmapping the file");
+ exit(EXIT_FAILURE);
+ }
+
+ map_abs[*fid] = map;
+
+ //printf("file length r: %d \n",filesize_abs[*fid]);
+ //printf("file map r: %d\n",*fid);
+
+}
+
+
+//void
+//FC_FUNC_(close_file_abs_map,CLOSE_FILE_ABS_MAP)(int * fid){
+void close_file_abs_map(int * fid){
+
+ /* Don't forget to free the mmapped memory
+ */
+ if (munmap(map_abs[*fid], filesize_abs[*fid]) == -1) {
+ perror("Error un-mmapping the file");
+ /* Decide here whether to close(fd) and exit() or not. Depends... */
+ }
+
+ /* Un-mmaping doesn't close the file, so we still need to do that.
+ */
+ close(map_fd_abs[*fid]);
+}
+
+
+//void
+//FC_FUNC_(write_abs_map,WRITE_ABS_MAP)(int *fid, char *buffer, int *length , int *index){
+void write_abs_map(int *fid, char *buffer, int *length , int *index){
+
+ char *map;
+ int offset;
+
+ map = map_abs[*fid];
+
+ // offset in bytes
+ offset = (*index -1 ) * (*length) ;
+
+ // copies buffer to map
+ memcpy( &map[offset], buffer ,*length );
+
+}
+
+//void
+//FC_FUNC_(read_abs_map,READ_ABS_MAP)(int *fid, char *buffer, int *length , int *index){
+void read_abs_map(int *fid, char *buffer, int *length , int *index){
+
+ char *map;
+ int offset;
+
+ map = map_abs[*fid];
+
+ // offset in bytes
+ offset = (*index -1 ) * (*length) ;
+
+ // copies map to buffer
+ memcpy( buffer, &map[offset], *length );
+
+}
+
+
+/*
+
+ wrapper functions
+
+ - for your preferred, optimized file i/o ;
+ e.g. uncomment // #define USE_MAP... in config.h to use mmap routines
+ or comment out (default) to use fopen/fwrite/fread functions
+
+ note: mmap functions should work fine for local harddisk directories, but can lead to
+ problems with global (e.g. NFS) directories
+
+ (on nehalem, Linux 2.6.18-164.11.1.el5 #1 SMP Wed Jan 20 10:04:55 EST 2010 x86_64 x86_64 x86_64 GNU/Linux
+ - mmap functions are about 20 % faster than conventional fortran, unformatted file i/o
+ - fwrite/fread function are about 12 % faster than conventional fortran, unformatted file i/o )
+
+ */
+
+void
+FC_FUNC_(open_file_abs_w,OPEN_FILE_ABS_W)(int *fid, char *filename,int *length, int *filesize) {
+
+#ifdef USE_MAP_FUNCTION
+ open_file_abs_w_map(fid,filename,length,filesize);
+#else
+ open_file_abs_w_fbin(fid,filename,length,filesize);
+#endif
+
+}
+
+void
+FC_FUNC_(open_file_abs_r,OPEN_FILE_ABS_R)(int *fid, char *filename,int *length, int *filesize) {
+
+#ifdef USE_MAP_FUNCTION
+ open_file_abs_r_map(fid,filename,length,filesize);
+#else
+ open_file_abs_r_fbin(fid,filename,length,filesize);
+#endif
+
+}
+
+void
+FC_FUNC_(close_file_abs,CLOSE_FILES_ABS)(int *fid) {
+
+#ifdef USE_MAP_FUNCTION
+ close_file_abs_map(fid);
+#else
+ close_file_abs_fbin(fid);
+#endif
+
+}
+
+void
+FC_FUNC_(write_abs,WRITE_ABS)(int *fid, char *buffer, int *length , int *index) {
+
+#ifdef USE_MAP_FUNCTION
+ write_abs_map(fid,buffer,length,index);
+#else
+ write_abs_fbin(fid,buffer,length,index);
+#endif
+
+}
+
+void
+FC_FUNC_(read_abs,READ_ABS)(int *fid, char *buffer, int *length , int *index) {
+
+#ifdef USE_MAP_FUNCTION
+ read_abs_map(fid,buffer,length,index);
+#else
+ read_abs_fbin(fid,buffer,length,index);
+#endif
+
+}
+
+
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in 2011-04-02 10:11:08 UTC (rev 18165)
@@ -67,6 +67,7 @@
$O/assemble_MPI_scalar.o \
$O/check_mesh_resolution.o \
$O/comp_source_time_function.o \
+ $O/compute_adj_source_frechet.o \
$O/compute_arrays_source.o \
$O/create_name_database.o \
$O/create_serial_name_database.o \
@@ -96,6 +97,7 @@
$O/sort_array_coordinates.o \
$O/utm_geo.o \
$O/write_VTK_data.o \
+ $O/write_c_binary.o \
$(EMPTY_MACRO)
# solver objects - no statically allocated arrays anymore
@@ -438,6 +440,9 @@
$O/define_derivation_matrices.o: $(SHARED)constants.h ${SHARED}/define_derivation_matrices.f90
${FCCOMPILE_CHECK} -c -o $O/define_derivation_matrices.o ${SHARED}/define_derivation_matrices.f90
+$O/compute_adj_source_frechet.o: $(SHARED)constants.h compute_adj_source_frechet.f90
+ ${FCCOMPILE_CHECK} -c -o $O/compute_adj_source_frechet.o compute_adj_source_frechet.f90
+
$O/compute_arrays_source.o: $(SHARED)constants.h ${SHARED}/compute_arrays_source.f90
${FCCOMPILE_CHECK} -c -o $O/compute_arrays_source.o ${SHARED}/compute_arrays_source.f90
Added: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_adj_source_frechet.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_adj_source_frechet.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_adj_source_frechet.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -0,0 +1,206 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 2 . 0
+! ---------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA and University of Pau / CNRS / INRIA
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! November 2010
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+
+! compute the integrated derivatives of source parameters (M_jk and X_s)
+
+ subroutine compute_adj_source_frechet(displ_s,Mxx,Myy,Mzz,Mxy,Mxz,Myz,eps_s,eps_m_s, &
+ hxir,hetar,hgammar,hpxir,hpetar,hpgammar, hprime_xx,hprime_yy,hprime_zz, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+
+ implicit none
+
+ include 'constants.h'
+
+ ! input
+ real(kind=CUSTOM_REAL) :: displ_s(NDIM,NGLLX,NGLLY,NGLLZ)
+ double precision :: Mxx, Myy, Mzz, Mxy, Mxz, Myz
+ ! output
+ real(kind=CUSTOM_REAL) :: eps_s(NDIM,NDIM), eps_m_s(NDIM)
+
+ ! auxilliary
+ double precision :: hxir(NGLLX), hetar(NGLLY), hgammar(NGLLZ), &
+ hpxir(NGLLX),hpetar(NGLLY),hpgammar(NGLLZ)
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+! local variables
+ real(kind=CUSTOM_REAL) :: tempx1l,tempx2l,tempx3l, tempy1l,tempy2l,tempy3l, &
+ tempz1l,tempz2l,tempz3l, hp1, hp2, hp3, &
+ xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl, &
+ duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl, &
+ xix_s,xiy_s,xiz_s,etax_s,etay_s,etaz_s,gammax_s,gammay_s,gammaz_s, &
+ hlagrange_xi, hlagrange_eta, hlagrange_gamma, hlagrange
+
+ real(kind=CUSTOM_REAL) :: eps(NDIM,NDIM), eps_array(NDIM,NDIM,NGLLX,NGLLY,NGLLZ), &
+ eps_m_array(NGLLX,NGLLY,NGLLZ)
+
+ integer i,j,k,l
+
+
+! first compute the strain at all the GLL points of the source element
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ tempx1l = 0._CUSTOM_REAL
+ tempx2l = 0._CUSTOM_REAL
+ tempx3l = 0._CUSTOM_REAL
+
+ tempy1l = 0._CUSTOM_REAL
+ tempy2l = 0._CUSTOM_REAL
+ tempy3l = 0._CUSTOM_REAL
+
+ tempz1l = 0._CUSTOM_REAL
+ tempz2l = 0._CUSTOM_REAL
+ tempz3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ tempx1l = tempx1l + displ_s(1,l,j,k)*hp1
+ tempy1l = tempy1l + displ_s(2,l,j,k)*hp1
+ tempz1l = tempz1l + displ_s(3,l,j,k)*hp1
+
+ hp2 = hprime_yy(j,l)
+ tempx2l = tempx2l + displ_s(1,i,l,k)*hp2
+ tempy2l = tempy2l + displ_s(2,i,l,k)*hp2
+ tempz2l = tempz2l + displ_s(3,i,l,k)*hp2
+
+ hp3 = hprime_zz(k,l)
+ tempx3l = tempx3l + displ_s(1,i,j,l)*hp3
+ tempy3l = tempy3l + displ_s(2,i,j,l)*hp3
+ tempz3l = tempz3l + displ_s(3,i,j,l)*hp3
+ enddo
+
+! dudx
+ xixl = xix(i,j,k)
+ xiyl = xiy(i,j,k)
+ xizl = xiz(i,j,k)
+ etaxl = etax(i,j,k)
+ etayl = etay(i,j,k)
+ etazl = etaz(i,j,k)
+ gammaxl = gammax(i,j,k)
+ gammayl = gammay(i,j,k)
+ gammazl = gammaz(i,j,k)
+
+ duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+! strain eps_jk
+ eps(1,1) = duxdxl
+ eps(1,2) = (duxdyl + duydxl) / 2
+ eps(1,3) = (duxdzl + duzdxl) / 2
+ eps(2,2) = duydyl
+ eps(2,3) = (duydzl + duzdyl) / 2
+ eps(3,3) = duzdzl
+ eps(2,1) = eps(1,2)
+ eps(3,1) = eps(1,3)
+ eps(3,2) = eps(2,3)
+
+ eps_array(:,:,i,j,k) = eps(:,:)
+
+! Mjk eps_jk
+ eps_m_array(i,j,k) = Mxx * eps(1,1) + Myy * eps(2,2) + Mzz * eps(3,3) + &
+ 2 * (Mxy * eps(1,2) + Mxz * eps(1,3) + Myz * eps(2,3))
+
+ enddo
+ enddo
+ enddo
+
+ ! interpolate the strain eps_s(:,:) from eps_array(:,:,i,j,k)
+ eps_s = 0.
+ xix_s = 0.; xiy_s = 0.; xiz_s = 0.
+ etax_s = 0.; etay_s = 0.; etaz_s = 0.
+ gammax_s = 0.; gammay_s = 0.; gammaz_s = 0.
+
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ hlagrange = hxir(i)*hetar(j)*hgammar(k)
+
+ eps_s(1,1) = eps_s(1,1) + eps_array(1,1,i,j,k)*hlagrange
+ eps_s(1,2) = eps_s(1,2) + eps_array(1,2,i,j,k)*hlagrange
+ eps_s(1,3) = eps_s(1,3) + eps_array(1,3,i,j,k)*hlagrange
+ eps_s(2,2) = eps_s(2,2) + eps_array(2,2,i,j,k)*hlagrange
+ eps_s(2,3) = eps_s(2,3) + eps_array(2,3,i,j,k)*hlagrange
+ eps_s(3,3) = eps_s(3,3) + eps_array(3,3,i,j,k)*hlagrange
+
+ xix_s = xix_s + xix(i,j,k)*hlagrange
+ xiy_s = xiy_s + xiy(i,j,k)*hlagrange
+ xiz_s = xiz_s + xiz(i,j,k)*hlagrange
+ etax_s = etax_s + etax(i,j,k)*hlagrange
+ etay_s = etay_s + etay(i,j,k)*hlagrange
+ etaz_s = etaz_s + etaz(i,j,k)*hlagrange
+ gammax_s = gammax_s + gammax(i,j,k)*hlagrange
+ gammay_s = gammay_s + gammay(i,j,k)*hlagrange
+ gammaz_s = gammaz_s + gammaz(i,j,k)*hlagrange
+
+ enddo
+ enddo
+ enddo
+
+! for completion purpose, not used in specfem3D.f90
+ eps_s(2,1) = eps_s(1,2)
+ eps_s(3,1) = eps_s(1,3)
+ eps_s(3,2) = eps_s(2,3)
+
+! compute the gradient of M_jk * eps_jk, and then interpolate it
+
+ eps_m_s = 0.
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ hlagrange_xi = hpxir(i)*hetar(j)*hgammar(k)
+ hlagrange_eta = hxir(i)*hpetar(j)*hgammar(k)
+ hlagrange_gamma = hxir(i)*hetar(j)*hpgammar(k)
+
+ eps_m_s(1) = eps_m_s(1) + eps_m_array(i,j,k) * (hlagrange_xi * xix_s &
+ + hlagrange_eta * etax_s + hlagrange_gamma * gammax_s)
+ eps_m_s(2) = eps_m_s(2) + eps_m_array(i,j,k) * (hlagrange_xi * xiy_s &
+ + hlagrange_eta * etay_s + hlagrange_gamma * gammay_s)
+ eps_m_s(3) = eps_m_s(3) + eps_m_array(i,j,k) * (hlagrange_xi * xiz_s &
+ + hlagrange_eta * etaz_s + hlagrange_gamma * gammaz_s)
+
+ enddo
+ enddo
+ enddo
+
+ end subroutine compute_adj_source_frechet
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -165,7 +165,7 @@
ibool,ispec_is_inner,phase_is_inner, &
abs_boundary_jacobian2Dw,abs_boundary_ijk,abs_boundary_ispec, &
num_abs_boundary_faces,rhostore,kappastore,ispec_is_acoustic, &
- SIMULATION_TYPE,SAVE_FORWARD,NSTEP,it,myrank,NGLOB_ADJOINT, &
+ SIMULATION_TYPE,SAVE_FORWARD,NSTEP,it,NGLOB_ADJOINT, &
b_potential_dot_dot_acoustic,b_reclen_potential, &
b_absorb_potential,b_num_abs_boundary_faces)
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -171,7 +171,7 @@
abs_boundary_ijk,abs_boundary_ispec, &
num_abs_boundary_faces, &
veloc,rho_vp,rho_vs, &
- ispec_is_elastic,SIMULATION_TYPE,myrank,SAVE_FORWARD, &
+ ispec_is_elastic,SIMULATION_TYPE,SAVE_FORWARD, &
NSTEP,it,NGLOB_ADJOINT,b_accel, &
b_num_abs_boundary_faces,b_reclen_field,b_absorb_field )
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_noDev.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_noDev.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -32,8 +32,10 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool,&
ATTENUATION,&
- one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,&
- NSPEC_ATTENUATION_AB,R_xx,R_yy,R_xy,R_xz,R_yz, &
+ one_minus_sum_beta,factor_common, &
+ alphaval,betaval,gammaval,&
+ NSPEC_ATTENUATION_AB, &
+ R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilondev_xx,epsilondev_yy,epsilondev_xy,&
epsilondev_xz,epsilondev_yz,epsilon_trace_over_3, &
ANISOTROPY,NSPEC_ANISO, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -31,7 +31,7 @@
ibool,ispec_is_inner,phase_is_inner, &
abs_boundary_jacobian2Dw,abs_boundary_ijk,abs_boundary_ispec, &
num_abs_boundary_faces,rhostore,kappastore,ispec_is_acoustic,&
- SIMULATION_TYPE,SAVE_FORWARD,NSTEP,it,myrank,NGLOB_ADJOINT, &
+ SIMULATION_TYPE,SAVE_FORWARD,NSTEP,it,NGLOB_ADJOINT, &
b_potential_dot_dot_acoustic,b_reclen_potential, &
b_absorb_potential,b_num_abs_boundary_faces)
@@ -61,7 +61,7 @@
! adjoint simulations
integer:: SIMULATION_TYPE
- integer:: NSTEP,it,myrank,NGLOB_ADJOINT
+ integer:: NSTEP,it,NGLOB_ADJOINT
integer:: b_num_abs_boundary_faces,b_reclen_potential
real(kind=CUSTOM_REAL),dimension(NGLOB_ADJOINT) :: b_potential_dot_dot_acoustic
real(kind=CUSTOM_REAL),dimension(NGLLSQUARE,b_num_abs_boundary_faces):: b_absorb_potential
@@ -70,15 +70,20 @@
! local parameters
real(kind=CUSTOM_REAL) :: rhol,cpl,jacobianw
integer :: ispec,iglob,i,j,k,iface,igll
- !adjoint locals
- integer:: reclen1,reclen2
+ !integer:: reclen1,reclen2
! adjoint simulations:
if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
+ ! reads in absorbing boundary
! the index NSTEP-it+1 is valid if b_displ is read in after the Newark scheme
- read(IOABS_AC,rec=NSTEP-it+1) reclen1,b_absorb_potential,reclen2
- if (reclen1 /= b_reclen_potential .or. reclen1 /= reclen2) &
- call exit_mpi(myrank,'Error reading absorbing contribution b_absorb_potential')
+
+ ! uses fortran routine
+ !read(IOABS_AC,rec=NSTEP-it+1) reclen1,b_absorb_potential,reclen2
+ !if (reclen1 /= b_reclen_potential .or. reclen1 /= reclen2) &
+ ! call exit_mpi(0,'Error reading absorbing contribution b_absorb_potential')
+ ! uses c routine for faster reading
+ call read_abs(1,b_absorb_potential,b_reclen_potential,NSTEP-it+1)
+
endif !adjoint
! absorbs absorbing-boundary surface using Sommerfeld condition (vanishing field in the outer-space)
@@ -128,7 +133,12 @@
enddo ! num_abs_boundary_faces
! adjoint simulations: stores absorbed wavefield part
- if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. num_abs_boundary_faces > 0 ) &
- write(IOABS_AC,rec=it) b_reclen_potential,b_absorb_potential,b_reclen_potential
+ if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. num_abs_boundary_faces > 0 ) then
+ ! writes out absorbing boundary value
+ ! uses fortran routine
+ !write(IOABS_AC,rec=it) b_reclen_potential,b_absorb_potential,b_reclen_potential
+ ! uses c routine
+ call write_abs(1,b_absorb_potential,b_reclen_potential,it)
+ endif
end subroutine compute_stacey_acoustic
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_elastic.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_elastic.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -34,7 +34,7 @@
abs_boundary_ijk,abs_boundary_ispec, &
num_abs_boundary_faces, &
veloc,rho_vp,rho_vs, &
- ispec_is_elastic,SIMULATION_TYPE,myrank,SAVE_FORWARD, &
+ ispec_is_elastic,SIMULATION_TYPE,SAVE_FORWARD, &
NSTEP,it,NGLOB_ADJOINT,b_accel, &
b_num_abs_boundary_faces,b_reclen_field,b_absorb_field)
@@ -67,7 +67,7 @@
! adjoint simulations
integer:: SIMULATION_TYPE
- integer:: NSTEP,it,myrank,NGLOB_ADJOINT
+ integer:: NSTEP,it,NGLOB_ADJOINT
integer:: b_num_abs_boundary_faces,b_reclen_field
real(kind=CUSTOM_REAL),dimension(NDIM,NGLLSQUARE,b_num_abs_boundary_faces):: b_absorb_field
@@ -77,17 +77,21 @@
! local parameters
real(kind=CUSTOM_REAL) vx,vy,vz,nx,ny,nz,tx,ty,tz,vn,jacobianw
integer :: ispec,iglob,i,j,k,iface,igll
+ !integer:: reclen1,reclen2
- !adjoint locals
- integer:: reclen1,reclen2
-
! adjoint simulations:
if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
+ ! reads in absorbing boundary
! the index NSTEP-it+1 is valid if b_displ is read in after the Newark scheme
- read(IOABS,rec=NSTEP-it+1) reclen1,b_absorb_field,reclen2
- if (reclen1 /= b_reclen_field .or. reclen1 /= reclen2) &
- call exit_mpi(myrank,'Error reading absorbing contribution b_absorb_field')
+
+ ! uses fortran routine
+ !read(IOABS,rec=NSTEP-it+1) reclen1,b_absorb_field,reclen2
+ !if (reclen1 /= b_reclen_field .or. reclen1 /= reclen2) &
+ ! call exit_mpi(0,'Error reading absorbing contribution b_absorb_field')
+ ! uses c routine for faster reading
+ call read_abs(0,b_absorb_field,b_reclen_field,NSTEP-it+1)
+
endif !adjoint
@@ -150,8 +154,13 @@
enddo
! adjoint simulations: stores absorbed wavefield part
- if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. num_abs_boundary_faces > 0 ) &
- write(IOABS,rec=it) b_reclen_field,b_absorb_field,b_reclen_field
-
+ if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. num_abs_boundary_faces > 0 ) then
+ ! writes out absorbing boundary value
+ ! uses fortran routine
+ !write(IOABS,rec=it) b_reclen_field,b_absorb_field,b_reclen_field
+ ! uses c routine
+ call write_abs(0,b_absorb_field,b_reclen_field,it)
+ endif
+
end subroutine compute_stacey_elastic
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -89,8 +89,8 @@
if( num_abs_boundary_faces > 0 .and. (SIMULATION_TYPE == 3 .or. &
(SIMULATION_TYPE == 1 .and. SAVE_FORWARD)) ) then
- if( ELASTIC_SIMULATION) close(IOABS)
- if( ACOUSTIC_SIMULATION) close(IOABS_AC)
+ if( ELASTIC_SIMULATION) call close_file_abs(0) ! close(IOABS)
+ if( ACOUSTIC_SIMULATION) call close_file_abs(1) ! close(IOABS_AC)
endif
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -279,13 +279,6 @@
NGLOB_ADJOINT = 1
endif
- ! strain/attenuation
- if( ATTENUATION .and. SIMULATION_TYPE == 3 ) then
- NSPEC_ATT_AND_KERNEL = NSPEC_AB
- else
- NSPEC_ATT_AND_KERNEL = 1
- endif
-
! moho boundary
if( SAVE_MOHO_MESH .and. SIMULATION_TYPE == 3 ) then
NSPEC_BOUN = NSPEC_AB
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -531,16 +531,30 @@
if (SIMULATION_TYPE == 3) then
! opens existing files
- open(unit=IOABS,file=trim(prname)//'absorb_field.bin',status='old',&
- action='read',form='unformatted',access='direct', &
- recl=b_reclen_field+2*4,iostat=ier )
- if( ier /= 0 ) call exit_mpi(myrank,'error opening proc***_absorb_field.bin file')
+
+ ! uses fortran routines for reading
+ !open(unit=IOABS,file=trim(prname)//'absorb_field.bin',status='old',&
+ ! action='read',form='unformatted',access='direct', &
+ ! recl=b_reclen_field+2*4,iostat=ier )
+ !if( ier /= 0 ) call exit_mpi(myrank,'error opening proc***_absorb_field.bin file')
+ ! uses c routines for faster reading
+ call open_file_abs_r(0,trim(prname)//'absorb_field.bin', &
+ len_trim(trim(prname)//'absorb_field.bin'), &
+ b_reclen_field*NSTEP)
+
else
! opens new file
- open(unit=IOABS,file=trim(prname)//'absorb_field.bin',status='unknown',&
- form='unformatted',access='direct',&
- recl=b_reclen_field+2*4,iostat=ier )
- if( ier /= 0 ) call exit_mpi(myrank,'error opening proc***_absorb_field.bin file')
+
+ ! uses fortran routines for writing
+ !open(unit=IOABS,file=trim(prname)//'absorb_field.bin',status='unknown',&
+ ! form='unformatted',access='direct',&
+ ! recl=b_reclen_field+2*4,iostat=ier )
+ !if( ier /= 0 ) call exit_mpi(myrank,'error opening proc***_absorb_field.bin file')
+ ! uses c routines for faster writing (file index 0 for acoutic domain file)
+ call open_file_abs_w(0,trim(prname)//'absorb_field.bin', &
+ len_trim(trim(prname)//'absorb_field.bin'), &
+ b_reclen_field*NSTEP)
+
endif
endif
@@ -554,16 +568,28 @@
if (SIMULATION_TYPE == 3) then
! opens existing files
- open(unit=IOABS_AC,file=trim(prname)//'absorb_potential.bin',status='old',&
- action='read',form='unformatted',access='direct', &
- recl=b_reclen_potential+2*4,iostat=ier )
- if( ier /= 0 ) call exit_mpi(myrank,'error opening proc***_absorb_potential.bin file')
+ ! uses fortran routines for reading
+ !open(unit=IOABS_AC,file=trim(prname)//'absorb_potential.bin',status='old',&
+ ! action='read',form='unformatted',access='direct', &
+ ! recl=b_reclen_potential+2*4,iostat=ier )
+ !if( ier /= 0 ) call exit_mpi(myrank,'error opening proc***_absorb_potential.bin file')
+ ! uses c routines for faster reading
+ call open_file_abs_r(1,trim(prname)//'absorb_potential.bin', &
+ len_trim(trim(prname)//'absorb_potential.bin'), &
+ b_reclen_potential*NSTEP)
+
else
! opens new file
- open(unit=IOABS_AC,file=trim(prname)//'absorb_potential.bin',status='unknown',&
- form='unformatted',access='direct',&
- recl=b_reclen_potential+2*4,iostat=ier )
- if( ier /= 0 ) call exit_mpi(myrank,'error opening proc***_absorb_potential.bin file')
+ ! uses fortran routines for writing
+ !open(unit=IOABS_AC,file=trim(prname)//'absorb_potential.bin',status='unknown',&
+ ! form='unformatted',access='direct',&
+ ! recl=b_reclen_potential+2*4,iostat=ier )
+ !if( ier /= 0 ) call exit_mpi(myrank,'error opening proc***_absorb_potential.bin file')
+ ! uses c routines for faster writing (file index 1 for acoutic domain file)
+ call open_file_abs_w(1,trim(prname)//'absorb_potential.bin', &
+ len_trim(trim(prname)//'absorb_potential.bin'), &
+ b_reclen_potential*NSTEP)
+
endif
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -518,14 +518,12 @@
if( ier /= 0 ) stop 'error allocating array b_request_send_vector_ext_mesh etc.'
! allocates attenuation solids
- if( ATTENUATION ) then
- allocate(b_R_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
- b_R_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
- b_R_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
- b_R_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
- b_R_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS),stat=ier)
- if( ier /= 0 ) stop 'error allocating array b_R_xx etc.'
- endif
+ allocate(b_R_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS), &
+ b_R_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS), &
+ b_R_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS), &
+ b_R_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS), &
+ b_R_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_R_xx etc.'
! note: these arrays are needed for attenuation and/or kernel computations
allocate(b_epsilondev_xx(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY), &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -735,14 +735,15 @@
double precision :: shape3D(NGNOD)
double precision :: xil,etal,gammal
double precision :: xmesh,ymesh,zmesh
-
real(kind=CUSTOM_REAL),dimension(NGNOD) :: xelm,yelm,zelm
+ integer :: ia,ispec,isource,irec,ier
+ character(len=256) :: filename,filename_new,system_command
- integer :: ia,ispec,isource,irec
-
if (myrank == 0) then
! vtk file
- open(IOVTK,file=trim(OUTPUT_FILES)//'/sr.vtk',status='unknown')
+ open(IOVTK,file=trim(OUTPUT_FILES)//'/sr.vtk',status='unknown',iostat=ier)
+ if( ier /= 0 ) stop 'error opening sr.vtk file'
+ ! vtk header
write(IOVTK,'(a)') '# vtk DataFile Version 2.0'
write(IOVTK,'(a)') 'Source and Receiver VTK file'
write(IOVTK,'(a)') 'ASCII'
@@ -859,6 +860,26 @@
if( myrank == 0 ) then
write(IOVTK,*)
close(IOVTK)
+
+ ! creates additional receiver and source files
+ ! extracts receiver locations
+ filename = trim(OUTPUT_FILES)//'/sr.vtk'
+ filename_new = trim(OUTPUT_FILES)//'/receiver.vtk'
+ write(system_command, &
+ "('awk ',a1,'{if(NR<5) print $0;if(NR==6)print ',a1,'POINTS',i6,' float',a1,';if(NR>5+',i6,')print $0}',a1,' < ',a,' > ',a)")&
+ "'",'"',nrec,'"',NSOURCES,"'",trim(filename),trim(filename_new)
+ call system(system_command)
+
+ ! extracts source locations
+ filename_new = trim(OUTPUT_FILES)//'/source.vtk'
+ write(system_command, &
+ "('awk ',a1,'{if(NR< 6 + ',i6,') print $0}END{print}',a1,' < ',a,' > ',a)")&
+ "'",NSOURCES,"'",trim(filename),trim(filename_new)
+ call system(system_command)
+
+
endif
+
+
end subroutine setup_sources_receivers_VTKfile
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -313,8 +313,6 @@
b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: b_epsilon_trace_over_3
- integer:: NSPEC_ATT_AND_KERNEL
-
! adjoint kernels
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rho_kl, mu_kl, kappa_kl, &
rhop_kl, beta_kl, alpha_kl
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2011-04-02 04:05:46 UTC (rev 18164)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2011-04-02 10:11:08 UTC (rev 18165)
@@ -81,9 +81,7 @@
real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable:: &
displ_element,veloc_element,accel_element
real(kind=CUSTOM_REAL),dimension(1):: dummy
- integer :: ipoin,ispec,iglob,ispec2D
- integer :: i,j,k,ier
- logical :: is_done
+ integer :: ipoin,ispec,iglob,ispec2D,ier
! allocate array for single elements
allocate( displ_element(NDIM,NGLLX,NGLLY,NGLLZ), &
@@ -166,38 +164,10 @@
! acoustic domains
if( ispec_is_acoustic(ispec) ) then
- ! velocity vector
- is_done = .false.
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- if( iglob == ibool(i,j,k,ispec) ) then
- ! norm of displacement
- store_val_ux_external_mesh(NGLLX*NGLLY*(ispec2D-1)+ipoin) = &
- max(store_val_ux_external_mesh(NGLLX*NGLLY*(ispec2D-1)+ipoin), &
- sqrt(displ_element(1,i,j,k)**2 &
- + displ_element(2,i,j,k)**2 &
- + displ_element(3,i,j,k)**2))
- ! norm of velocity
- store_val_uy_external_mesh(NGLLX*NGLLY*(ispec2D-1)+ipoin) = &
- max(store_val_uy_external_mesh(NGLLX*NGLLY*(ispec2D-1)+ipoin), &
- sqrt(veloc_element(1,i,j,k)**2 &
- + veloc_element(2,i,j,k)**2 &
- + veloc_element(3,i,j,k)**2))
- ! norm of acceleration
- store_val_uz_external_mesh(NGLLX*NGLLY*(ispec2D-1)+ipoin) = &
- max(store_val_uz_external_mesh(NGLLX*NGLLY*(ispec2D-1)+ipoin), &
- sqrt(accel_element(1,i,j,k)**2 &
- + accel_element(2,i,j,k)**2 &
- + accel_element(3,i,j,k)**2))
- is_done = .true.
- exit
- endif
- enddo
- if( is_done ) exit
- enddo
- if( is_done ) exit
- enddo
+ ! sets velocity vector with maximum norm of wavefield values
+ call wmo_get_max_vector(ispec,ispec2D,ipoin, &
+ displ_element,veloc_element,accel_element, &
+ NGLLX*NGLLY)
endif
enddo
@@ -223,38 +193,10 @@
! acoustic domains
if( ispec_is_acoustic(ispec) ) then
- ! velocity vector
- is_done = .false.
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- if( iglob == ibool(i,j,k,ispec) ) then
- ! norm of displacement
- store_val_ux_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = &
- max(store_val_ux_external_mesh(NGNOD2D*(ispec2D-1)+ipoin), &
- sqrt(displ_element(1,i,j,k)**2 &
- + displ_element(2,i,j,k)**2 &
- + displ_element(3,i,j,k)**2))
- ! norm of velocity
- store_val_uy_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = &
- max(store_val_uy_external_mesh(NGNOD2D*(ispec2D-1)+ipoin), &
- sqrt(veloc_element(1,i,j,k)**2 &
- + veloc_element(2,i,j,k)**2 &
- + veloc_element(3,i,j,k)**2))
- ! norm of acceleration
- store_val_uz_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = &
- max(store_val_uz_external_mesh(NGNOD2D*(ispec2D-1)+ipoin), &
- sqrt(accel_element(1,i,j,k)**2 &
- + accel_element(2,i,j,k)**2 &
- + accel_element(3,i,j,k)**2))
- is_done = .true.
- exit
- endif
- enddo
- if( is_done ) exit
- enddo
- if( is_done ) exit
- enddo
+ ! sets velocity vector with maximum norm of wavefield values
+ call wmo_get_max_vector(ispec,ispec2D,ipoin, &
+ displ_element,veloc_element,accel_element, &
+ NGNOD2D)
endif
enddo
endif
@@ -322,7 +264,59 @@
end subroutine wmo_create_shakemap_em
+!================================================================
+ subroutine wmo_get_max_vector(ispec,ispec2D,ipoin, &
+ displ_element,veloc_element,accel_element, &
+ narraydim)
+
+ ! put into this separate routine to make compilation faster
+
+ use specfem_par,only: NDIM,ibool
+ use specfem_par_movie
+ implicit none
+
+ integer :: ispec,ispec2D,ipoin,narraydim
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: &
+ displ_element,veloc_element,accel_element
+
+ ! local parameters
+ integer :: i,j,k,iglob
+ logical :: is_done
+
+ is_done = .false.
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ if( iglob == ibool(i,j,k,ispec) ) then
+ ! norm of displacement
+ store_val_ux_external_mesh(narraydim*(ispec2D-1)+ipoin) = &
+ max(store_val_ux_external_mesh(narraydim*(ispec2D-1)+ipoin), &
+ sqrt(displ_element(1,i,j,k)**2 &
+ + displ_element(2,i,j,k)**2 &
+ + displ_element(3,i,j,k)**2))
+ ! norm of velocity
+ store_val_uy_external_mesh(narraydim*(ispec2D-1)+ipoin) = &
+ max(store_val_uy_external_mesh(narraydim*(ispec2D-1)+ipoin), &
+ sqrt(veloc_element(1,i,j,k)**2 &
+ + veloc_element(2,i,j,k)**2 &
+ + veloc_element(3,i,j,k)**2))
+ ! norm of acceleration
+ store_val_uz_external_mesh(narraydim*(ispec2D-1)+ipoin) = &
+ max(store_val_uz_external_mesh(narraydim*(ispec2D-1)+ipoin), &
+ sqrt(accel_element(1,i,j,k)**2 &
+ + accel_element(2,i,j,k)**2 &
+ + accel_element(3,i,j,k)**2))
+ ! not really needed, but could be used to check...
+ is_done = .true.
+ return
+ endif
+ enddo
+ enddo
+ enddo
+
+ end subroutine wmo_get_max_vector
+
!================================================================
subroutine wmo_create_movie_surface_em()
@@ -337,8 +331,7 @@
real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable:: veloc_element
real(kind=CUSTOM_REAL),dimension(1):: dummy
- integer :: ispec2D,ispec,ipoin,iglob,i,j,k,ier
- logical :: is_done
+ integer :: ispec2D,ispec,ipoin,iglob,ier
! allocate array for single elements
allocate( veloc_element(NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
@@ -393,23 +386,10 @@
! acoustic pressure potential
if( ispec_is_acoustic(ispec) ) then
- ! velocity vector
- is_done = .false.
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- if( iglob == ibool(i,j,k,ispec) ) then
- store_val_ux_external_mesh(NGLLX*NGLLY*(ispec2D-1)+ipoin) = veloc_element(1,i,j,k)
- store_val_uy_external_mesh(NGLLX*NGLLY*(ispec2D-1)+ipoin) = veloc_element(2,i,j,k)
- store_val_uz_external_mesh(NGLLX*NGLLY*(ispec2D-1)+ipoin) = veloc_element(3,i,j,k)
- is_done = .true.
- exit
- endif
- enddo
- if( is_done ) exit
- enddo
- if( is_done ) exit
- enddo
+ ! puts velocity values into storage array
+ call wmo_get_vel_vector(ispec,ispec2D,ipoin, &
+ veloc_element, &
+ NGLLX*NGLLY)
endif
enddo
else
@@ -425,23 +405,10 @@
! acoustic pressure potential
if( ispec_is_acoustic(ispec) ) then
- ! velocity vector
- is_done = .false.
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- if( iglob == ibool(i,j,k,ispec) ) then
- store_val_ux_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = veloc_element(1,i,j,k)
- store_val_uy_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = veloc_element(2,i,j,k)
- store_val_uz_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = veloc_element(3,i,j,k)
- is_done = .true.
- exit
- endif
- enddo
- if( is_done ) exit
- enddo
- if( is_done ) exit
- enddo
+ ! puts velocity values into storage array
+ call wmo_get_vel_vector(ispec,ispec2D,ipoin, &
+ veloc_element, &
+ NGNOD2D)
endif
enddo
endif
@@ -517,7 +484,45 @@
end subroutine wmo_create_movie_surface_em
+!================================================================
+ subroutine wmo_get_vel_vector(ispec,ispec2D,ipoin, &
+ veloc_element, &
+ narraydim)
+
+ ! put into this separate routine to make compilation faster
+
+ use specfem_par,only: NDIM,ibool
+ use specfem_par_movie
+ implicit none
+
+ integer :: ispec,ispec2D,ipoin,narraydim
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: &
+ veloc_element
+
+ ! local parameters
+ integer :: i,j,k,iglob
+ logical :: is_done
+
+ ! velocity vector
+ is_done = .false.
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ if( iglob == ibool(i,j,k,ispec) ) then
+ store_val_ux_external_mesh(narraydim*(ispec2D-1)+ipoin) = veloc_element(1,i,j,k)
+ store_val_uy_external_mesh(narraydim*(ispec2D-1)+ipoin) = veloc_element(2,i,j,k)
+ store_val_uz_external_mesh(narraydim*(ispec2D-1)+ipoin) = veloc_element(3,i,j,k)
+ is_done = .true.
+ return
+ endif
+ enddo
+ enddo
+ enddo
+
+ end subroutine wmo_get_vel_vector
+
+
!=====================================================================
subroutine wmo_movie_surface_output_o()
@@ -534,7 +539,6 @@
real(kind=CUSTOM_REAL),dimension(1) :: dummy
integer :: ispec,ipoin,iglob,i,j,k,ier
integer :: imin,imax,jmin,jmax,kmin,kmax,iface,igll,iloc
- logical :: is_done
! allocate array for single elements
allocate( val_element(NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
@@ -632,23 +636,8 @@
! acoustic pressure potential
if( ispec_is_acoustic(ispec) ) then
- ! velocity vector
- is_done = .false.
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- if( iglob == ibool(i,j,k,ispec) ) then
- store_val_ux_external_mesh(ipoin) = val_element(1,i,j,k)
- store_val_uy_external_mesh(ipoin) = val_element(2,i,j,k)
- store_val_uz_external_mesh(ipoin) = val_element(3,i,j,k)
- is_done = .true.
- exit
- endif
- enddo
- if( is_done ) exit
- enddo
- if( is_done ) exit
- enddo
+ ! stores values from element
+ call wmo_get_val_elem(ispec,ipoin,val_element)
endif
enddo
@@ -685,23 +674,8 @@
! acoustic pressure potential
if( ispec_is_acoustic(ispec) ) then
- ! velocity vector
- is_done = .false.
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- if( iglob == ibool(i,j,k,ispec) ) then
- store_val_ux_external_mesh(ipoin) = val_element(1,i,j,k)
- store_val_uy_external_mesh(ipoin) = val_element(2,i,j,k)
- store_val_uz_external_mesh(ipoin) = val_element(3,i,j,k)
- is_done = .true.
- exit
- endif
- enddo
- if( is_done ) exit
- enddo
- if( is_done ) exit
- enddo
+ ! stores values from element
+ call wmo_get_val_elem(ispec,ipoin,val_element)
endif
enddo ! iloc
@@ -775,7 +749,42 @@
end subroutine wmo_movie_surface_output_o
+!================================================================
+ subroutine wmo_get_val_elem(ispec,ipoin,val_element)
+
+ ! put into this separate routine to make compilation faster
+
+ use specfem_par,only: NDIM,ibool
+ use specfem_par_movie
+ implicit none
+
+ integer :: ispec,ipoin
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: &
+ val_element
+
+ ! local parameters
+ integer :: i,j,k,iglob
+ logical :: is_done
+
+ ! velocity vector
+ is_done = .false.
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ if( iglob == ibool(i,j,k,ispec) ) then
+ store_val_ux_external_mesh(ipoin) = val_element(1,i,j,k)
+ store_val_uy_external_mesh(ipoin) = val_element(2,i,j,k)
+ store_val_uz_external_mesh(ipoin) = val_element(3,i,j,k)
+ is_done = .true.
+ return
+ endif
+ enddo
+ enddo
+ enddo
+
+ end subroutine wmo_get_val_elem
+
!=====================================================================
subroutine wmo_create_shakemap_o()
@@ -794,7 +803,6 @@
integer :: ipoin,ispec,iglob
integer :: imin,imax,jmin,jmax,kmin,kmax,iface,igll,iloc
integer :: i,j,k,ier
- logical :: is_done
! allocate array for single elements
allocate( displ_element(NDIM,NGLLX,NGLLY,NGLLZ), &
@@ -828,7 +836,6 @@
ibool,rhostore)
endif
-
! save all points for high resolution, or only four corners for low resolution
if(USE_HIGHRES_FOR_MOVIES) then
do igll = 1, NGLLSQUARE
@@ -855,30 +862,8 @@
! acoustic domains
if( ispec_is_acoustic(ispec) ) then
- ! velocity vector
- is_done = .false.
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- if( iglob == ibool(i,j,k,ispec) ) then
- ! horizontal displacement
- store_val_ux_external_mesh(ipoin) = max(store_val_ux_external_mesh(ipoin),&
- abs(displ_element(1,i,j,k)),abs(displ_element(2,i,j,k)))
- ! horizontal velocity
- store_val_uy_external_mesh(ipoin) = max(store_val_uy_external_mesh(ipoin),&
- abs(veloc_element(1,i,j,k)),abs(veloc_element(2,i,j,k)))
- ! horizontal acceleration
- store_val_uz_external_mesh(ipoin) = max(store_val_uz_external_mesh(ipoin),&
- abs(accel_element(1,i,j,k)),abs(accel_element(2,i,j,k)))
-
- is_done = .true.
- exit
- endif
- enddo
- if( is_done ) exit
- enddo
- if( is_done ) exit
- enddo
+ ! stores maximum values
+ call wmo_get_max_vector_o(ispec,ipoin,displ_element,veloc_element,accel_element)
endif
enddo
@@ -915,30 +900,8 @@
! acoustic domains
if( ispec_is_acoustic(ispec) ) then
- ! velocity vector
- is_done = .false.
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- if( iglob == ibool(i,j,k,ispec) ) then
- ! horizontal displacement
- store_val_ux_external_mesh(ipoin) = max(store_val_ux_external_mesh(ipoin),&
- abs(displ_element(1,i,j,k)),abs(displ_element(2,i,j,k)))
- ! horizontal velocity
- store_val_uy_external_mesh(ipoin) = max(store_val_uy_external_mesh(ipoin),&
- abs(veloc_element(1,i,j,k)),abs(veloc_element(2,i,j,k)))
- ! horizontal acceleration
- store_val_uz_external_mesh(ipoin) = max(store_val_uz_external_mesh(ipoin),&
- abs(accel_element(1,i,j,k)),abs(accel_element(2,i,j,k)))
-
- is_done = .true.
- exit
- endif
- enddo
- if( is_done ) exit
- enddo
- if( is_done ) exit
- enddo
+ ! stores maximum values
+ call wmo_get_max_vector_o(ispec,ipoin,displ_element,veloc_element,accel_element)
endif
enddo
@@ -1008,6 +971,49 @@
end subroutine wmo_create_shakemap_o
+!================================================================
+
+ subroutine wmo_get_max_vector_o(ispec,ipoin,displ_element,veloc_element,accel_element)
+
+ ! put into this separate routine to make compilation faster
+
+ use specfem_par,only: NDIM,ibool
+ use specfem_par_movie
+ implicit none
+
+ integer :: ispec,ipoin
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: &
+ displ_element,veloc_element,accel_element
+
+ ! local parameters
+ integer :: i,j,k,iglob
+ logical :: is_done
+
+ ! velocity vector
+ is_done = .false.
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ if( iglob == ibool(i,j,k,ispec) ) then
+ ! horizontal displacement
+ store_val_ux_external_mesh(ipoin) = max(store_val_ux_external_mesh(ipoin),&
+ abs(displ_element(1,i,j,k)),abs(displ_element(2,i,j,k)))
+ ! horizontal velocity
+ store_val_uy_external_mesh(ipoin) = max(store_val_uy_external_mesh(ipoin),&
+ abs(veloc_element(1,i,j,k)),abs(veloc_element(2,i,j,k)))
+ ! horizontal acceleration
+ store_val_uz_external_mesh(ipoin) = max(store_val_uz_external_mesh(ipoin),&
+ abs(accel_element(1,i,j,k)),abs(accel_element(2,i,j,k)))
+
+ is_done = .true.
+ return
+ endif
+ enddo
+ enddo
+ enddo
+
+ end subroutine wmo_get_max_vector_o
+
!=====================================================================
subroutine wmo_movie_volume_output()
More information about the CIG-COMMITS
mailing list