[cig-commits] r20688 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src

[surendra at geodynamics.org]

Author: surendra
Date: 2012-09-05 07:49:33 -0700 (Wed, 05 Sep 2012)
New Revision: 20688

Modified:
   seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/check_mesh_resolution.f90
Log:
Fixed a bug in checking highest frequency resolved based on development version of SPECFEM3D

Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/check_mesh_resolution.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/check_mesh_resolution.f90	2012-09-05 09:27:25 UTC (rev 20687)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/check_mesh_resolution.f90	2012-09-05 14:49:33 UTC (rev 20688)
@@ -63,6 +63,24 @@
   real(kind=CUSTOM_REAL),parameter :: NELEM_PER_WAVELENGTH = 1.5
   logical :: has_vs_zero
 
+!Surendra
+  ! corners indices of reference cube faces
+  ! shapes of arrays below
+  integer, parameter :: NGNOD_EIGHT_CORNERS = 8
+  integer,dimension(2),parameter :: corner_shape = (/3,NGNOD_EIGHT_CORNERS/)
+  integer,dimension(3,NGNOD_EIGHT_CORNERS),parameter :: corner_ijk = &
+    reshape((/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ, &
+      NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),corner_shape)
+  ! local parameters
+  real(kind=CUSTOM_REAL) :: x0,y0,z0
+  integer :: icorner,jcorner
+  real(kind=CUSTOM_REAL) :: elemsize_min,elemsize_max,elemsize_min_glob,elemsize_max_glob
+  real(kind=CUSTOM_REAL) :: avg_distance
+  ! number of points per minimum wavelength for minimum period estimate
+  real(kind=CUSTOM_REAL),parameter :: NPTS_PER_WAVELENGTH = 5
+
+
+
 ! initializations
   if( DT <= 0.0d0) then
     DT_PRESENT = .false.
@@ -167,6 +185,45 @@
     distance_min_glob = min( distance_min_glob, distance_min)
     distance_max_glob = max( distance_max_glob, distance_max)
     
+  ! initializes
+  elemsize_min = HUGEVAL
+  elemsize_max = -HUGEVAL
+
+  ! loops over corners
+  do icorner=1,NGNOD_EIGHT_CORNERS
+    i = corner_ijk(1,icorner)
+    j = corner_ijk(2,icorner)
+    k = corner_ijk(3,icorner)
+
+    iglob_a = ibool(i,j,k,ispec)
+    x0 = xstore(iglob_a)
+    y0 = ystore(iglob_a)
+    z0 = zstore(iglob_a)
+
+    ! loops over all other corners
+    do jcorner = icorner+1,NGNOD_EIGHT_CORNERS
+      i = corner_ijk(1,jcorner)
+      j = corner_ijk(2,jcorner)
+      k = corner_ijk(3,jcorner)
+
+      ! coordinates
+      iglob_b = ibool(i,j,k,ispec)
+
+      ! distances between points
+      if( iglob_a /= iglob_b) then
+        dx = sqrt( ( x0 - xstore(iglob_b) )**2 &
+                  + ( y0 - ystore(iglob_b) )**2 &
+                  + ( z0 - zstore(iglob_b) )**2 )
+        if( dx < elemsize_min) elemsize_min = dx
+        if( dx > elemsize_max) elemsize_max = dx
+      endif
+
+    enddo
+  enddo
+
+    elemsize_min_glob = min( elemsize_min_glob, elemsize_min)
+    elemsize_max_glob = max( elemsize_max_glob, elemsize_max)
+
     ! courant number
     if( DT_PRESENT ) then
       cmax = max( vpmax,vsmax ) * DT / distance_min    
@@ -177,8 +234,26 @@
     dt_suggested = COURANT_SUGGESTED * distance_min / max( vpmax,vsmax )
     dt_suggested_glob = min( dt_suggested_glob, dt_suggested)          
                             
+   ! estimation of minimum period resolved
+    ! based on average GLL distance within element and minimum velocity
+    !
+    ! rule of thumb (Komatitsch et al. 2005):
+    ! "average number of points per minimum wavelength in an element should be around 5."
+
+    ! average distance between GLL points within this element
+    avg_distance = elemsize_max / ( NGLLX - 1 )  ! since NGLLX = NGLLY = NGLLZ
+
+    ! largest possible minimum period such that number of points per minimum wavelength
+    ! npts = ( min(vpmin,vsmin)  * pmax ) / avg_distance  is about ~ NPTS_PER_WAVELENGTH
+    !
+    ! note: obviously, this estimation depends on the choice of points per wavelength
+    !          which is empirical at the moment.
+    !          also, keep in mind that the minimum period is just an estimation and
+    !          there is no such sharp cut-off period for valid synthetics.
+    !          seismograms become just more and more inaccurate for periods shorter than this estimate.
+    pmax = avg_distance / min( vpmin,vsmin ) * NPTS_PER_WAVELENGTH
     ! estimation of minimum period resolved
-    pmax = distance_max / min( vpmin,vsmin ) * NELEM_PER_WAVELENGTH  
+    !pmax = distance_max / min( vpmin,vsmin ) * NELEM_PER_WAVELENGTH  
     pmax_glob = max(pmax_glob,pmax)
     
   enddo
@@ -297,4 +372,4 @@
   
   
   end subroutine
-  
\ No newline at end of file
+