[cig-commits] r19334 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Jan 8 16:19:07 PST 2012
Author: dkomati1
Date: 2012-01-08 16:19:07 -0800 (Sun, 08 Jan 2012)
New Revision: 19334
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/checkgrid.F90
Log:
modified the wavelength sampling histogram to compute separately the S wave histogram in solid regions and the P wave histogram in fluid regions
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/checkgrid.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/checkgrid.F90 2012-01-07 20:46:21 UTC (rev 19333)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/checkgrid.F90 2012-01-09 00:19:07 UTC (rev 19334)
@@ -136,7 +136,7 @@
double precision :: distance_min_glob,distance_max_glob
double precision :: courant_stability_max_glob,lambdaPImin_glob,lambdaPImax_glob,&
lambdaPIImin_glob,lambdaPIImax_glob,lambdaSmin_glob,lambdaSmax_glob, &
- lambdaSmin_glob_with_P_in_fluid,lambdaSmax_glob_with_P_in_fluid
+ lambdaPmin_in_fluid_histo_glob,lambdaPmax_in_fluid_histo_glob
double precision :: xmin_glob, xmax_glob, zmin_glob, zmax_glob
#endif
@@ -155,12 +155,13 @@
integer :: is,ir,in,nnum
! for histogram of number of points per wavelength
- double precision :: lambdaSmin_P_in_fluid_histo,lambdaSmax_P_in_fluid_histo,lambdaPmin_histo,lambdaPmax_histo, &
+ double precision :: lambdaPmin_in_fluid_histo,lambdaPmax_in_fluid_histo,lambdaSmin_histo,lambdaSmax_histo, &
min_nb_of_points_per_wavelength,max_nb_of_points_per_wavelength,nb_of_points_per_wavelength, &
scaling_factor,scaling_factor_S,scaling_factor_P
+ logical :: any_fluid_histo,any_fluid_histo_glob
integer, parameter :: NCLASSES = 20
integer, dimension(0:NCLASSES-1) :: classes_wavelength,classes_wavelength_all
- integer :: iclass,nspec_all,ipass
+ integer :: iclass,nspec_all,ipass,nspec_counted,nspec_counted_all,nspec_counted_all_solid,nspec_counted_all_fluid
logical :: create_wavelength_histogram
double precision :: current_percent,total_percent
@@ -240,9 +241,11 @@
lambdaPIImax = 0
endif
- lambdaSmin_P_in_fluid_histo = HUGEVAL
- lambdaSmax_P_in_fluid_histo = -HUGEVAL
+ lambdaPmin_in_fluid_histo = HUGEVAL
+ lambdaPmax_in_fluid_histo = -HUGEVAL
+ any_fluid_histo = .false.
+
do ispec=1,nspec
material = kmato(ispec)
@@ -312,11 +315,11 @@
vpImax = max(vpImax,cpIloc)
! ignore acoustic and elastic regions with cpII = 0
- if(cpIIloc > 0.0001d0) vpIImin = min(vpIImin,cpIIloc)
+ if(cpIIloc > 1.d-20) vpIImin = min(vpIImin,cpIIloc)
vpIImax = max(vpIImax,cpIIloc)
! ignore fluid regions with Vs = 0
- if(csloc > 0.0001d0) vsmin = min(vsmin,csloc)
+ if(csloc > 1.d-20) vsmin = min(vsmin,csloc)
vsmax = max(vsmax,csloc)
densmin = min(densmin,denst)
@@ -353,29 +356,26 @@
courant_stability_number_max = max(courant_stability_number_max, &
vpImax_local * deltat / (distance_min_local * percent_GLL(NGLLX)))
-! ignore fluid regions with Vs = 0
- if(vsmin_local > 0.0001d0) then
+! check if fluid region with Vs = 0
+ if(vsmin_local > 1.d-20) then
lambdaSmin = min(lambdaSmin,vsmin_local / (distance_max_local / (NGLLX - 1)))
lambdaSmax = max(lambdaSmax,vsmin_local / (distance_max_local / (NGLLX - 1)))
- lambdaSmin_P_in_fluid_histo = min(lambdaSmin_P_in_fluid_histo,vsmin_local / (distance_max_local / (NGLLX - 1)))
- lambdaSmax_P_in_fluid_histo = max(lambdaSmax_P_in_fluid_histo,vsmin_local / (distance_max_local / (NGLLX - 1)))
else
- lambdaSmin_P_in_fluid_histo = min(lambdaSmin_P_in_fluid_histo,vpImin_local / (distance_max_local / (NGLLX - 1)))
- lambdaSmax_P_in_fluid_histo = max(lambdaSmax_P_in_fluid_histo,vpImin_local / (distance_max_local / (NGLLX - 1)))
+ any_fluid_histo = .true.
+ lambdaPmin_in_fluid_histo = min(lambdaPmin_in_fluid_histo,vpImin_local / (distance_max_local / (NGLLX - 1)))
+ lambdaPmax_in_fluid_histo = max(lambdaPmax_in_fluid_histo,vpImin_local / (distance_max_local / (NGLLX - 1)))
endif
lambdaPImin = min(lambdaPImin,vpImin_local / (distance_max_local / (NGLLX - 1)))
lambdaPImax = max(lambdaPImax,vpImin_local / (distance_max_local / (NGLLX - 1)))
- if(cpIIloc > 0.0001d0) then
+ if(cpIIloc > 1.d-20) then
lambdaPIImin = min(lambdaPIImin,vpIImin_local / (distance_max_local / (NGLLX - 1)))
lambdaPIImax = max(lambdaPIImax,vpIImin_local / (distance_max_local / (NGLLX - 1)))
endif
enddo
- any_elastic_glob = any_elastic
- any_poroelastic_glob = any_poroelastic
#ifdef USE_MPI
call MPI_ALLREDUCE (vpImin, vpImin_glob, 1, MPI_DOUBLE_PRECISION, &
MPI_MIN, MPI_COMM_WORLD, ier)
@@ -411,14 +411,16 @@
MPI_MIN, MPI_COMM_WORLD, ier)
call MPI_ALLREDUCE (lambdaSmax, lambdaSmax_glob, 1, MPI_DOUBLE_PRECISION, &
MPI_MAX, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (lambdaSmin_P_in_fluid_histo, lambdaSmin_glob_with_P_in_fluid, 1, MPI_DOUBLE_PRECISION, &
+ call MPI_ALLREDUCE (lambdaPmin_in_fluid_histo, lambdaPmin_in_fluid_histo_glob, 1, MPI_DOUBLE_PRECISION, &
MPI_MIN, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (lambdaSmax_P_in_fluid_histo, lambdaSmax_glob_with_P_in_fluid, 1, MPI_DOUBLE_PRECISION, &
+ call MPI_ALLREDUCE (lambdaPmax_in_fluid_histo, lambdaPmax_in_fluid_histo_glob, 1, MPI_DOUBLE_PRECISION, &
MPI_MAX, MPI_COMM_WORLD, ier)
call MPI_ALLREDUCE (any_elastic, any_elastic_glob, 1, MPI_LOGICAL, &
MPI_LOR, MPI_COMM_WORLD, ier)
call MPI_ALLREDUCE (any_poroelastic, any_poroelastic_glob, 1, MPI_LOGICAL, &
MPI_LOR, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (any_fluid_histo, any_fluid_histo_glob, 1, MPI_LOGICAL, &
+ MPI_LOR, MPI_COMM_WORLD, ier)
vpImin = vpImin_glob
vpImax = vpImax_glob
vpIImin = vpIImin_glob
@@ -436,8 +438,12 @@
lambdaPIImax = lambdaPIImax_glob
lambdaSmin = lambdaSmin_glob
lambdaSmax = lambdaSmax_glob
- lambdaSmin_P_in_fluid_histo = lambdaSmin_glob_with_P_in_fluid
- lambdaSmax_P_in_fluid_histo = lambdaSmax_glob_with_P_in_fluid
+ lambdaPmin_in_fluid_histo = lambdaPmin_in_fluid_histo_glob
+ lambdaPmax_in_fluid_histo = lambdaPmax_in_fluid_histo_glob
+#else
+ any_elastic_glob = any_elastic
+ any_poroelastic_glob = any_poroelastic
+ any_fluid_histo_glob = any_fluid_histo
#endif
if ( myrank == 0 ) then
@@ -492,15 +498,11 @@
write(IOUT,*) '----'
! for histogram
- lambdaSmin_P_in_fluid_histo = lambdaSmin_P_in_fluid_histo/(2.5d0*f0max)
- lambdaSmax_P_in_fluid_histo = lambdaSmax_P_in_fluid_histo/(2.5d0*f0max)
+ lambdaPmin_in_fluid_histo = lambdaPmin_in_fluid_histo/(2.5d0*f0max)
+ lambdaPmax_in_fluid_histo = lambdaPmax_in_fluid_histo/(2.5d0*f0max)
- if(lambdaPIImin <= ZERO) then
- lambdaPmin_histo = lambdaPImin/(2.5d0*f0max)
- else
- lambdaPmin_histo = min(lambdaPImin,lambdaPIImin)/(2.5d0*f0max)
- endif
- lambdaPmax_histo = max(lambdaPImax,lambdaPIImax)/(2.5d0*f0max)
+ lambdaSmin_histo = lambdaSmin/(2.5d0*f0max)
+ lambdaSmax_histo = lambdaSmax/(2.5d0*f0max)
create_wavelength_histogram = .true.
@@ -512,10 +514,10 @@
endif
#ifdef USE_MPI
- call MPI_BCAST(lambdaSmin_P_in_fluid_histo,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
- call MPI_BCAST(lambdaSmax_P_in_fluid_histo,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
- call MPI_BCAST(lambdaPmin_histo,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
- call MPI_BCAST(lambdaPmax_histo,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(lambdaPmin_in_fluid_histo,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(lambdaPmax_in_fluid_histo,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(lambdaSmin_histo,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(lambdaSmax_histo,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
call MPI_BCAST(f0max,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
call MPI_BCAST(create_wavelength_histogram,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ier)
#endif
@@ -529,15 +531,24 @@
! create statistics about mesh sampling (number of points per wavelength)
-! first pass is for S wave sampling, second pass is for P wave sampling
+ nspec_counted_all_solid = 0
+ nspec_counted_all_fluid = 0
+
+! first pass is for S wave sampling in solid, second pass is for P wave sampling in fluid
do ipass = 1,2
+ nspec_counted = 0
+
if(ipass == 1) then
- min_nb_of_points_per_wavelength = lambdaSmin_P_in_fluid_histo
- max_nb_of_points_per_wavelength = lambdaSmax_P_in_fluid_histo
+ min_nb_of_points_per_wavelength = lambdaSmin_histo
+ max_nb_of_points_per_wavelength = lambdaSmax_histo
+! do not create this histogram if the model is entirely fluid
+ if(.not. any_elastic_glob .and. .not. any_poroelastic_glob) cycle
else
- min_nb_of_points_per_wavelength = lambdaPmin_histo
- max_nb_of_points_per_wavelength = lambdaPmax_histo
+! do not create this histogram if the model is entirely solid
+ if(.not. any_fluid_histo_glob) cycle
+ min_nb_of_points_per_wavelength = lambdaPmin_in_fluid_histo
+ max_nb_of_points_per_wavelength = lambdaPmax_in_fluid_histo
endif
! when the grid is regular and the medium is homogeneous, the minimum and the maximum are equal
@@ -637,32 +648,47 @@
if(ipass == 1) then
-! ignore fluid regions with Vs = 0
- if(vsmin_local > 0.0001d0) then
+! in first pass, only solid regions, thus ignore fluid regions with Vs = 0
+ if(vsmin_local > 1.d-20) then
nb_of_points_per_wavelength = vsmin_local / (distance_max_local / (NGLLX - 1))
- else
- nb_of_points_per_wavelength = vpImin_local / (distance_max_local / (NGLLX - 1))
+
+ nspec_counted = nspec_counted + 1
+
+ nb_of_points_per_wavelength = nb_of_points_per_wavelength/(2.5d0*f0max)
+
+! store number of points per wavelength in histogram
+ iclass = int((nb_of_points_per_wavelength - min_nb_of_points_per_wavelength) / &
+ (max_nb_of_points_per_wavelength - min_nb_of_points_per_wavelength) * dble(NCLASSES))
+ if(iclass < 0) iclass = 0
+ if(iclass > NCLASSES-1) iclass = NCLASSES-1
+ classes_wavelength(iclass) = classes_wavelength(iclass) + 1
+
endif
else
- if(vpIImin_local <= ZERO) then
- nb_of_points_per_wavelength = vpImin_local / (distance_max_local / (NGLLX - 1))
- else
- nb_of_points_per_wavelength = min(vpImin_local,vpIImin_local) / (distance_max_local / (NGLLX - 1))
- endif
+! in second pass, only fluid regions, thus ignore solid regions with Vs > 0
+ if(abs(vsmin_local) < 1.d-20) then
+ if(vpIImin_local <= ZERO) then
+ nb_of_points_per_wavelength = vpImin_local / (distance_max_local / (NGLLX - 1))
+ else
+ nb_of_points_per_wavelength = min(vpImin_local,vpIImin_local) / (distance_max_local / (NGLLX - 1))
+ endif
- endif
+ nspec_counted = nspec_counted + 1
- nb_of_points_per_wavelength = nb_of_points_per_wavelength/(2.5d0*f0max)
+ nb_of_points_per_wavelength = nb_of_points_per_wavelength/(2.5d0*f0max)
! store number of points per wavelength in histogram
- iclass = int((nb_of_points_per_wavelength - min_nb_of_points_per_wavelength) / &
- (max_nb_of_points_per_wavelength - min_nb_of_points_per_wavelength) * dble(NCLASSES))
- if(iclass < 0) iclass = 0
- if(iclass > NCLASSES-1) iclass = NCLASSES-1
- classes_wavelength(iclass) = classes_wavelength(iclass) + 1
+ iclass = int((nb_of_points_per_wavelength - min_nb_of_points_per_wavelength) / &
+ (max_nb_of_points_per_wavelength - min_nb_of_points_per_wavelength) * dble(NCLASSES))
+ if(iclass < 0) iclass = 0
+ if(iclass > NCLASSES-1) iclass = NCLASSES-1
+ classes_wavelength(iclass) = classes_wavelength(iclass) + 1
+ endif
+ endif
+
enddo
#ifdef USE_MPI
@@ -673,38 +699,53 @@
#ifdef USE_MPI
call MPI_REDUCE(nspec, nspec_all, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ier)
+ call MPI_REDUCE(nspec_counted, nspec_counted_all, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ier)
#else
+ nspec_counted_all = nspec_counted
nspec_all = nspec
#endif
+ if(ipass == 1) then
+ nspec_counted_all_solid = nspec_counted_all
+ else
+ nspec_counted_all_fluid = nspec_counted_all
+ endif
+
! create histogram of wavelength and save in Gnuplot file
if (myrank == 0) then
write(IOUT,*)
+ write(IOUT,*) '-----------------------------------------'
+ write(IOUT,*)
if(ipass == 1) then
- write(IOUT,*) 'histogram of min number of points per S wavelength'
- write(IOUT,*) ' (P wavelength in acoustic regions)'
+ write(IOUT,*) 'histogram of min number of points per S wavelength in solid regions:'
+ write(IOUT,*)
+ write(IOUT,*) 'there are ',nspec_counted_all,' elements out of ',nspec_all,' in solid regions'
else
- write(IOUT,*) 'histogram of min number of points per P wavelength'
+ write(IOUT,*) 'histogram of min number of points per P wavelength in fluid regions:'
+ write(IOUT,*)
+ write(IOUT,*) 'there are ',nspec_counted_all,' elements out of ',nspec_all,' in fluid regions'
endif
- write(IOUT,*) '(too small: poor resolution of calculations -'
+ write(IOUT,*) ' (i.e., ',sngl(100.d0*nspec_counted_all/dble(nspec_all)),'% of the total)'
+ write(IOUT,*)
+ write(IOUT,*) '(too small = poor resolution of calculations -'
write(IOUT,*) ' too big = wasting memory and CPU time)'
- write(IOUT,*) '(threshold value is around 4.5 points per wavelength'
- write(IOUT,*) ' in elastic media and 5.5 in acoustic media):'
+ write(IOUT,*) '(threshold value is around 4.5 points per S wavelength'
+ write(IOUT,*) ' in elastic regions and 5.5 per P wavelength in fluid regions):'
write(IOUT,*)
total_percent = 0.
scaling_factor = max_nb_of_points_per_wavelength - min_nb_of_points_per_wavelength
if(ipass == 1) then
- open(unit=14,file='points_per_wavelength_histogram_S.txt',status='unknown')
+ open(unit=14,file='points_per_wavelength_histogram_S_in_solid.txt',status='unknown')
scaling_factor_S = scaling_factor
else
- open(unit=14,file='points_per_wavelength_histogram_P.txt',status='unknown')
+ open(unit=14,file='points_per_wavelength_histogram_P_in_fluid.txt',status='unknown')
scaling_factor_P = scaling_factor
endif
do iclass = 0,NCLASSES-1
- current_percent = 100.*dble(classes_wavelength_all(iclass))/dble(nspec_all)
+ current_percent = 100.*dble(classes_wavelength_all(iclass))/dble(nspec_counted_all)
total_percent = total_percent + current_percent
write(IOUT,*) sngl(min_nb_of_points_per_wavelength + scaling_factor*iclass/dble(NCLASSES)),' - ', &
sngl(min_nb_of_points_per_wavelength + scaling_factor*(iclass+1)/dble(NCLASSES)),classes_wavelength_all(iclass), &
@@ -717,6 +758,9 @@
if(total_percent < 99.9d0 .or. total_percent > 100.1d0) then
write(IOUT,*) 'total percentage = ',total_percent,' %'
stop 'total percentage should be 100%'
+ else
+ write(IOUT,*)
+ write(IOUT,*) 'total percentage = ',total_percent,' %'
endif
endif ! of if myrank == 0
@@ -729,29 +773,34 @@
open(unit=14,file='plot_points_per_wavelength_histogram.gnu',status='unknown')
write(14,*) 'set term x11'
- write(14,*) '#set term gif'
- write(14,*) '#set output "points_per_wavelength_histogram_S.gif"'
- write(14,*)
- write(14,*) 'set boxwidth ',real(scaling_factor_S/NCLASSES)
- write(14,*) 'set xlabel "Range of min number of points per S wavelength"'
- write(14,*) 'set ylabel "Percentage of elements (%)"'
- write(14,*) 'plot "points_per_wavelength_histogram_S.txt" with boxes'
- write(14,*) 'pause -1 "hit any key..."'
+ if(nspec_counted_all_solid > 0) then
+ write(14,*) '#set term gif'
+ write(14,*) '#set output "points_per_wavelength_histogram_S_in_solid.gif"'
+ write(14,*)
+ write(14,*) 'set boxwidth ',real(scaling_factor_S/NCLASSES)
+ write(14,*) 'set xlabel "Range of min number of points per S wavelength in solid"'
+ write(14,*) 'set ylabel "Percentage of elements (%)"'
+ write(14,*) 'plot "points_per_wavelength_histogram_S_in_solid.txt" with boxes'
+ write(14,*) 'pause -1 "hit any key..."'
+ endif
- write(14,*) '#set term gif'
- write(14,*) '#set output "points_per_wavelength_histogram_P.gif"'
- write(14,*)
- write(14,*) 'set boxwidth ',real(scaling_factor_P/NCLASSES)
- write(14,*) 'set xlabel "Range of min number of points per P wavelength"'
- write(14,*) 'set ylabel "Percentage of elements (%)"'
- write(14,*) 'plot "points_per_wavelength_histogram_P.txt" with boxes'
- write(14,*) 'pause -1 "hit any key..."'
+ if(nspec_counted_all_fluid > 0) then
+ write(14,*) '#set term gif'
+ write(14,*) '#set output "points_per_wavelength_histogram_P_in_fluid.gif"'
+ write(14,*)
+ write(14,*) 'set boxwidth ',real(scaling_factor_P/NCLASSES)
+ write(14,*) 'set xlabel "Range of min number of points per P wavelength in fluid"'
+ write(14,*) 'set ylabel "Percentage of elements (%)"'
+ write(14,*) 'plot "points_per_wavelength_histogram_P_in_fluid.txt" with boxes'
+ write(14,*) 'pause -1 "hit any key..."'
+ endif
close(14)
write(IOUT,*)
- write(IOUT,*) 'total number of elements = ',nspec_all
write(IOUT,*)
+ write(IOUT,*) 'total number of elements in fluid and solid regions = ',nspec_all
+ write(IOUT,*)
endif ! of if myrank == 0
@@ -1159,7 +1208,7 @@
!
!---- open PostScript file
!
- if(any_elastic_glob .or. any_poroelastic) then
+ if(any_elastic_glob .or. any_poroelastic_glob) then
open(unit=24,file='OUTPUT_FILES/mesh_S_wave_dispersion.ps',status='unknown')
else
open(unit=24,file='OUTPUT_FILES/mesh_P_wave_dispersion.ps',status='unknown')
@@ -1441,7 +1490,7 @@
if(any_elastic_glob .or. any_poroelastic_glob) then
! ignore fluid regions with Vs = 0
- if(vsmin_local > 0.0001d0) then
+ if(vsmin_local > 1.d-20) then
lambdaS_local = vsmin_local / (distance_max_local / (NGLLX - 1))
More information about the CIG-COMMITS
mailing list