[cig-commits] [commit] devel: - added a print statement to prepare_color_image if the code detects that the image is too big - increased the image size limit to 10000 pixels in each direction - added a small Fortran program to create a pseudo plane wave source based on tapering with a Hamming window - fixed the "Number of absorbing elements" print statement bug: this should not be printed when PML is on (84cb404)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Feb 13 08:15:14 PST 2014
Repository : ssh://geoshell/specfem2d
On branch : devel
Link : https://github.com/geodynamics/specfem2d/compare/594a714d2be70c6eedd3f07b2dcbb97364c4e08c...3bae725501ad5989ffbc8b45021f2471cb801bad
>---------------------------------------------------------------
commit 84cb404ce0de3d3928e5f38294ea3531aab2b2de
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Thu Feb 13 17:12:51 2014 +0100
- added a print statement to prepare_color_image if the code detects that the image is too big
- increased the image size limit to 10000 pixels in each direction
- added a small Fortran program to create a pseudo plane wave source based on tapering with a Hamming window
- fixed the "Number of absorbing elements" print statement bug: this should not be printed when PML is on
>---------------------------------------------------------------
84cb404ce0de3d3928e5f38294ea3531aab2b2de
...ate_pressure_plane_wave_with_Hamming_window.f90 | 73 ++++++++++++++++++++++
setup/constants.h.in | 2 +-
src/specfem2D/prepare_color_image.F90 | 10 ++-
src/specfem2D/read_databases.F90 | 11 ++--
src/specfem2D/specfem2D.F90 | 2 +-
5 files changed, 88 insertions(+), 10 deletions(-)
diff --git a/UTILS/create_pressure_plane_wave_with_Hamming_window.f90 b/UTILS/create_pressure_plane_wave_with_Hamming_window.f90
new file mode 100644
index 0000000..09e5d85
--- /dev/null
+++ b/UTILS/create_pressure_plane_wave_with_Hamming_window.f90
@@ -0,0 +1,73 @@
+
+ program create_Hamming_window
+
+! create a "pseudo plane wave" source file using tapering based on a Hamming window
+
+! Dimitri Komatitsch, CNRS Marseille, France, February 2014,
+! based on discussions with Paul Cristini, also from CNRS Marseille, France.
+
+ implicit none
+
+! pi
+ double precision, parameter :: PI = 3.141592653589793d0
+
+ integer, parameter :: NSOURCES = 1000
+
+! the plane wave will extend in the vertical direction from zmin to zmax
+! and its amplitude in between will vary as a Hamming apodization window
+ double precision, parameter :: xs_center = -0.1d0
+ double precision, parameter :: zs_center = 0.07d0
+ double precision, parameter :: zs_size = 0.43d0
+ double precision, parameter :: zs_min = zs_center - zs_size/2.d0
+ double precision, parameter :: zs_max = zs_center + zs_size/2.d0
+
+! angle in degrees by which we rotate the plane wave
+ double precision, parameter :: angle_rotate_source = -45.d0 * (PI / 180.d0)
+
+ double precision, parameter :: factor_max = 1.d10
+
+ integer :: isource
+
+ double precision :: x,hamming,xs,zs,xval,zval,xprime,zprime
+
+ do isource = 1,NSOURCES
+
+! Hamming apodization window
+! see e.g. http://docs.scipy.org/doc/numpy/reference/generated/numpy.hamming.html
+! and http://www.mathworks.fr/fr/help/signal/ref/hamming.html
+ x = dble(isource - 1) / dble(NSOURCES - 1)
+ hamming = 0.54d0 - 0.46d0*cos(2*PI*x)
+
+ xs = xs_center
+ zs = zs_min + x * zs_size
+
+! subtract the position of the center of rotation
+ xval = xs - xs_center
+ zval = zs - zs_center
+
+! rotate it
+ xprime = xval*cos(angle_rotate_source) - zval*sin(angle_rotate_source)
+ zprime = xval*sin(angle_rotate_source) + zval*cos(angle_rotate_source)
+
+! add the position of the center of rotation back
+ xprime = xprime + xs_center
+ zprime = zprime + zs_center
+
+ write(*,*) '# source ',isource
+ write(*,*) 'source_surf = .false.'
+ write(*,*) 'xs = ',xprime
+ write(*,*) 'zs = ',zprime
+ write(*,*) 'source_type = 1'
+ write(*,*) 'time_function_type = 1'
+ write(*,*) 'f0 = 1.d6'
+ write(*,*) 't0 = 0.d0'
+ write(*,*) 'angleforce = 0.0'
+ write(*,*) 'Mxx = 1.'
+ write(*,*) 'Mzz = 1.'
+ write(*,*) 'Mxz = 0.'
+ write(*,*) 'factor = ',factor_max * hamming
+
+ enddo
+
+ end program create_Hamming_window
+
diff --git a/setup/constants.h.in b/setup/constants.h.in
index 90684fd..b03d6b5 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -31,7 +31,7 @@
integer, parameter :: NGLLZ = NGLLX
! maximum size of the X and Z directions of a JPEG image generated by the display routine
- integer, parameter :: NX_NZ_IMAGE_MAX = 2500
+ integer, parameter :: NX_NZ_IMAGE_MAX = 10000
! further reduce cache misses inner/outer in two passes in the case of an MPI simulation
! this flag is ignored in the case of a serial simulation
diff --git a/src/specfem2D/prepare_color_image.F90 b/src/specfem2D/prepare_color_image.F90
index 4b21441..be1f697 100644
--- a/src/specfem2D/prepare_color_image.F90
+++ b/src/specfem2D/prepare_color_image.F90
@@ -124,10 +124,16 @@
if (NZ_IMAGE_color > 65534) &
call exit_MPI('output image too big: NZ_IMAGE_color > 65534; increase factor_subsample_image in DATA/Par_file.')
- if (NX_IMAGE_color > NX_NZ_IMAGE_MAX) call exit_MPI( &
+ if (NX_IMAGE_color > NX_NZ_IMAGE_MAX) then
+ print *,'NX_IMAGE_color,NX_NZ_IMAGE_MAX = ',NX_IMAGE_color,NX_NZ_IMAGE_MAX
+ call exit_MPI( &
'output image too big: NX_IMAGE_color > NX_NZ_IMAGE_MAX; increase factor_subsample_image or change NX_NZ_IMAGE_MAX.')
- if (NZ_IMAGE_color > NX_NZ_IMAGE_MAX) call exit_MPI( &
+ endif
+ if (NZ_IMAGE_color > NX_NZ_IMAGE_MAX) then
+ print *,'NZ_IMAGE_color,NX_NZ_IMAGE_MAX = ',NZ_IMAGE_color,NX_NZ_IMAGE_MAX
+ call exit_MPI( &
'output image too big: NZ_IMAGE_color > NX_NZ_IMAGE_MAX; increase factor_subsample_image or change NX_NZ_IMAGE_MAX.')
+ endif
end subroutine prepare_color_image_init
diff --git a/src/specfem2D/read_databases.F90 b/src/specfem2D/read_databases.F90
index cf73772..4b28f37 100644
--- a/src/specfem2D/read_databases.F90
+++ b/src/specfem2D/read_databases.F90
@@ -518,7 +518,7 @@
!---- print element group main parameters
if (myrank == 0 .and. ipass == 1) then
write(IOUT,107)
- write(IOUT,207) nspec,ngnod,NGLLX,NGLLZ,NGLLX*NGLLZ,pointsdisp,numat,nelemabs
+ write(IOUT,207) nspec,ngnod,NGLLX,NGLLZ,NGLLX*NGLLZ,pointsdisp,numat
endif
! output formats
@@ -530,8 +530,7 @@
'Number of points in Y-direction . . . (NGLLZ) =',i7,/5x, &
'Number of points per element. . .(NGLLX*NGLLZ) =',i7,/5x, &
'Number of points for display . . .(pointsdisp) =',i7,/5x, &
- 'Number of element material sets . . . (numat) =',i7,/5x, &
- 'Number of absorbing elements . . . .(nelemabs) =',i7)
+ 'Number of element material sets . . . (numat) =',i7)
end subroutine read_databases_coorg_elem
@@ -683,7 +682,7 @@
ibegin_edge3,iend_edge3,ibegin_edge4,iend_edge4, &
numabs,codeabs,typeabs,perm,antecedent_list, &
nspec_left,nspec_right,nspec_bottom,nspec_top, &
- ib_right,ib_left,ib_bottom,ib_top)
+ ib_right,ib_left,ib_bottom,ib_top,PML_BOUNDARY_CONDITIONS)
! reads in absorbing edges
@@ -700,7 +699,7 @@
ibegin_edge3,iend_edge3,ibegin_edge4,iend_edge4,ibegin_edge2,iend_edge2
logical, dimension(4,nelemabs) :: codeabs
integer, dimension(nelemabs) :: typeabs
- logical :: anyabs
+ logical :: anyabs,PML_BOUNDARY_CONDITIONS
integer, dimension(nspec) :: perm,antecedent_list
integer :: nspec_left,nspec_right,nspec_bottom,nspec_top
@@ -840,7 +839,7 @@
endif
#endif
- if (myrank == 0 .and. ipass == 1) then
+ if (myrank == 0 .and. ipass == 1 .and. .not. PML_BOUNDARY_CONDITIONS) then
write(IOUT,*)
write(IOUT,*) 'Number of absorbing elements: ',nelemabs_tot
write(IOUT,*) ' nspec_left = ',nspec_left_tot
diff --git a/src/specfem2D/specfem2D.F90 b/src/specfem2D/specfem2D.F90
index c4e7112..74d7bf8 100644
--- a/src/specfem2D/specfem2D.F90
+++ b/src/specfem2D/specfem2D.F90
@@ -1550,7 +1550,7 @@
ibegin_edge3,iend_edge3,ibegin_edge4,iend_edge4, &
numabs,codeabs,typeabs,perm,antecedent_list, &
nspec_left,nspec_right,nspec_bottom,nspec_top, &
- ib_right,ib_left,ib_bottom,ib_top)
+ ib_right,ib_left,ib_bottom,ib_top,PML_BOUNDARY_CONDITIONS)
if(anyabs .and. (.not. PML_BOUNDARY_CONDITIONS))then
STACEY_BOUNDARY_CONDITIONS = .true.
More information about the CIG-COMMITS
mailing list