[cig-commits] r22417 - in seismo/3D/SPECFEM3D/trunk/src: generate_databases shared
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Tue Jun 25 05:15:59 PDT 2013
Author: danielpeter
Date: 2013-06-25 05:15:59 -0700 (Tue, 25 Jun 2013)
New Revision: 22417
Modified:
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_absorbing_boundary.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/setup_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90
Log:
bug fix for simulations without absorbing boundaries; slight modifications and user output splitup in check_mesh_resolution.f90
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2013-06-25 12:01:17 UTC (rev 22416)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2013-06-25 12:15:59 UTC (rev 22417)
@@ -257,11 +257,24 @@
max_interface_size_ext_mesh,ibool_interfaces_ext_mesh, &
SAVE_MESH_FILES,ANISOTROPY)
-! call fault_save_arrays_test(prname) ! for debugging
- call fault_save_arrays(prname)
+! saves faults
+ if( ANY_FAULT ) then
+ call sync_all()
+ if( myrank == 0) then
+ write(IMAIN,*) ' ...saving fault databases'
+ call flush_IMAIN()
+ endif
+ ! call fault_save_arrays_test(prname) ! for debugging
+ call fault_save_arrays(prname)
+ endif
! saves moho surface
if( SAVE_MOHO_MESH ) then
+ call sync_all()
+ if( myrank == 0) then
+ write(IMAIN,*) ' ...saving Moho surfaces'
+ call flush_IMAIN()
+ endif
call crm_save_moho()
endif
@@ -273,6 +286,10 @@
! checks the mesh, stability and resolved period
call sync_all()
+ if( myrank == 0) then
+ write(IMAIN,*) ' ...checking mesh resolution'
+ call flush_IMAIN()
+ endif
if( POROELASTIC_SIMULATION ) then
!chris: check for poro: At the moment cpI & cpII are for eta=0
@@ -291,6 +308,11 @@
! saves binary mesh files for attenuation
if( ATTENUATION ) then
+ call sync_all()
+ if( myrank == 0) then
+ write(IMAIN,*) ' ...saving attenuation databases'
+ call flush_IMAIN()
+ endif
call get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENUATION_RATIO, &
mustore,rho_vs,kappastore,rho_vp,qkappa_attenuation_store,qmu_attenuation_store, &
ispec_is_elastic,min_resolved_period,prname,FULL_ATTENUATION_SOLID)
@@ -309,35 +331,38 @@
deallocate(xstore_dummy,ystore_dummy,zstore_dummy)
endif
+ if( ACOUSTIC_SIMULATION) then
+ deallocate(rmass_acoustic)
+ deallocate(rmass_acoustic_interface)
+ endif
+
+ if( ELASTIC_SIMULATION ) then
+ deallocate(rmass)
+ endif
+
+ if( POROELASTIC_SIMULATION) then
+ deallocate(rmass_solid_poroelastic,rmass_fluid_poroelastic)
+ endif
+
if(STACEY_ABSORBING_CONDITIONS)then
if( ELASTIC_SIMULATION ) then
- deallocate(rmass)
deallocate(rmassx,rmassy,rmassz)
endif
if( ACOUSTIC_SIMULATION) then
- deallocate(rmass_acoustic)
- deallocate(rmass_acoustic_interface)
deallocate(rmassz_acoustic)
endif
endif
if(PML_CONDITIONS)then
if( ELASTIC_SIMULATION ) then
- deallocate(rmass)
if(PML_CONDITIONS)then
if(ACOUSTIC_SIMULATION)then
write(IOUT)rmass_elastic_interface
endif
endif
endif
- if( ACOUSTIC_SIMULATION) then
- deallocate(rmass_acoustic)
- endif
endif
- if( POROELASTIC_SIMULATION) then
- deallocate(rmass_solid_poroelastic,rmass_fluid_poroelastic)
- endif
end subroutine create_regions_mesh
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_absorbing_boundary.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_absorbing_boundary.f90 2013-06-25 12:01:17 UTC (rev 22416)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_absorbing_boundary.f90 2013-06-25 12:15:59 UTC (rev 22417)
@@ -493,9 +493,8 @@
enddo
enddo
endif
- endif
- if(PML_CONDITIONS)then
+ else if(PML_CONDITIONS)then
if( .not. PML_INSTEAD_OF_FREE_SURFACE ) then
! stores free surface
! sets face infos
@@ -546,7 +545,25 @@
enddo
enddo
endif
+
+ else
+ ! stores free surface
+ ! sets face infos
+ ifree = ifree + 1
+ free_surface_ispec(ifree) = ispec
+
+ ! gll points -- assuming NGLLX = NGLLY = NGLLZ
+ igllfree = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ igllfree = igllfree+1
+ free_surface_ijk(:,igllfree,ifree) = ijk_face(:,i,j)
+ free_surface_jacobian2Dw(igllfree,ifree) = jacobian2Dw_face(i,j)
+ free_surface_normal(:,igllfree,ifree) = normal_face(:,i,j)
+ enddo
+ enddo
endif
+
enddo
! checks counters
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2013-06-25 12:01:17 UTC (rev 22416)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2013-06-25 12:15:59 UTC (rev 22417)
@@ -384,6 +384,7 @@
!
subroutine save_arrays_solver_files(nspec,nglob,ibool)
+ use generate_databases_par, only: myrank
use create_regions_mesh_ext_par
implicit none
@@ -400,8 +401,13 @@
integer :: j,inum
character(len=256) :: filename
- logical,parameter :: DEBUG = .true.
+ logical,parameter :: DEBUG = .false.
+ if( myrank == 0) then
+ write(IMAIN,*) ' saving mesh files for AVS, OpenDX, Paraview'
+ call flush_IMAIN()
+ endif
+
! mesh arrays used for example in combine_vol_data.f90
!--- x coordinate
open(unit=27,file=prname(1:len_trim(prname))//'x.bin',status='unknown',form='unformatted',iostat=ier)
@@ -492,6 +498,12 @@
! VTK file output
if( DEBUG ) then
+ call sync_all()
+ if( myrank == 0) then
+ write(IMAIN,*) ' saving debugging mesh files'
+ call flush_IMAIN()
+ endif
+
! saves free surface points
if( num_free_surface_faces > 0 ) then
! saves free surface interface points
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/setup_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/setup_mesh.f90 2013-06-25 12:01:17 UTC (rev 22416)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/setup_mesh.f90 2013-06-25 12:15:59 UTC (rev 22417)
@@ -84,9 +84,10 @@
call max_all_dp(max_elevation,max_elevation_all)
if(myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) 'min and max of topography included in mesh in m is ',min_elevation_all,' ',max_elevation_all
- write(IMAIN,*)
+ write(IMAIN,*)
+ write(IMAIN,*) 'min and max of topography included in mesh in m is ',min_elevation_all,' ',max_elevation_all
+ write(IMAIN,*)
+ call flush_IMAIN()
endif
! clean-up
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2013-06-25 12:01:17 UTC (rev 22416)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2013-06-25 12:15:59 UTC (rev 22417)
@@ -62,7 +62,7 @@
integer :: myrank
integer :: NSPEC_AB_global_min,NSPEC_AB_global_max,NSPEC_AB_global_sum
integer :: NGLOB_AB_global_min,NGLOB_AB_global_max,NGLOB_AB_global_sum
- integer :: ispec,sizeprocs
+ integer :: ispec !,sizeprocs
!********************************************************************************
@@ -126,6 +126,44 @@
tmp2(:) = 0.0
endif
+!! DK DK May 2009: added this to print the minimum and maximum number of elements
+!! DK DK May 2009: and points in the CUBIT + SCOTCH mesh
+ call min_all_i(NSPEC_AB,NSPEC_AB_global_min)
+ call max_all_i(NSPEC_AB,NSPEC_AB_global_max)
+ call sum_all_i(NSPEC_AB,NSPEC_AB_global_sum)
+
+ call min_all_i(NGLOB_AB,NGLOB_AB_global_min)
+ call max_all_i(NGLOB_AB,NGLOB_AB_global_max)
+ call sum_all_i(NGLOB_AB,NGLOB_AB_global_sum)
+
+! outputs infos
+ if ( myrank == 0 ) then
+ write(IMAIN,*)
+ write(IMAIN,*) '********'
+ write(IMAIN,*) 'minimum and maximum number of elements'
+ write(IMAIN,*) 'and points in the CUBIT + SCOTCH mesh:'
+ write(IMAIN,*)
+ write(IMAIN,*) 'NSPEC_global_min = ',NSPEC_AB_global_min
+ write(IMAIN,*) 'NSPEC_global_max = ',NSPEC_AB_global_max
+ write(IMAIN,*) 'NSPEC_global_max / NSPEC_global_min imbalance = ',sngl(dble(NSPEC_AB_global_max) / dble(NSPEC_AB_global_min)),&
+ ' = ',sngl((dble(NSPEC_AB_global_max) / dble(NSPEC_AB_global_min) - 1.d0) * 100.d0),' %'
+ write(IMAIN,*) 'NSPEC_global_sum = ',NSPEC_AB_global_sum
+ write(IMAIN,*)
+ write(IMAIN,*) 'NGLOB_global_min = ',NGLOB_AB_global_min
+ write(IMAIN,*) 'NGLOB_global_max = ',NGLOB_AB_global_max
+ write(IMAIN,*) 'NGLOB_global_max / NGLOB_global_min imbalance = ',sngl(dble(NGLOB_AB_global_max) / dble(NGLOB_AB_global_min)),&
+ ' = ',sngl((dble(NGLOB_AB_global_max) / dble(NGLOB_AB_global_min) - 1.d0) * 100.d0),' %'
+ write(IMAIN,*) 'NGLOB_global_sum = ',NGLOB_AB_global_sum
+ write(IMAIN,*)
+ write(IMAIN,*) 'If you have elements of a single type (all acoustic, all elastic, all poroelastic, and without CPML)'
+ write(IMAIN,*) 'in the whole mesh, then there should be no significant imbalance in the above numbers.'
+ write(IMAIN,*) 'Otherwise, it is normal to have imbalance in elements and points because the domain decomposer'
+ write(IMAIN,*) 'compensates for the different cost of different elements by partitioning them unevenly among processes.'
+ write(IMAIN,*) '********'
+ write(IMAIN,*)
+ call flush_IMAIN()
+ endif
+
! checks courant number & minimum resolved period for each grid cell
do ispec=1,NSPEC_AB
@@ -208,21 +246,6 @@
enddo
-! determines global min/max values from all cpu partitions
- if( DT_PRESENT ) then
- ! courant number
- cmax = cmax_glob
- call max_all_cr(cmax,cmax_glob)
- endif
-
- ! minimum period
- pmax = pmax_glob
- call max_all_cr(pmax,pmax_glob)
-
- ! time step
- dt_suggested = dt_suggested_glob
- call min_all_cr(dt_suggested,dt_suggested_glob)
-
! Vp velocity
vpmin = vpmin_glob
vpmax = vpmax_glob
@@ -238,18 +261,31 @@
call max_all_cr(vsmax,vsmax_glob)
! checks velocities
- if( vpmin_glob <= 0.0_CUSTOM_REAL ) then
- call exit_mpi(myrank,"error: vp minimum velocity")
+ if( myrank == 0 ) then
+ if( vpmin_glob <= 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: vp minimum velocity")
+ endif
+ if( vpmax_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: vp maximum velocity")
+ endif
+ if( vsmin_glob < 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: vs minimum velocity")
+ endif
+ if( vsmax_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: vs maximum velocity")
+ endif
endif
- if( vpmax_glob >= HUGEVAL ) then
- call exit_mpi(myrank,"error: vp maximum velocity")
+
+! outputs infos
+ if ( myrank == 0 ) then
+ write(IMAIN,*)
+ write(IMAIN,*) '********'
+ write(IMAIN,*) 'Model: P velocity min,max = ',vpmin_glob,vpmax_glob
+ write(IMAIN,*) 'Model: S velocity min,max = ',vsmin_glob,vsmax_glob
+ write(IMAIN,*) '********'
+ write(IMAIN,*)
+ call flush_IMAIN()
endif
- if( vsmin_glob < 0.0_CUSTOM_REAL ) then
- call exit_mpi(myrank,"error: vs minimum velocity")
- endif
- if( vsmax_glob >= HUGEVAL ) then
- call exit_mpi(myrank,"error: vs maximum velocity")
- endif
! GLL point distance
distance_min = distance_min_glob
@@ -257,6 +293,28 @@
call min_all_cr(distance_min,distance_min_glob)
call max_all_cr(distance_max,distance_max_glob)
+ ! element size
+ elemsize_min = elemsize_min_glob
+ elemsize_max = elemsize_max_glob
+ call min_all_cr(elemsize_min,elemsize_min_glob)
+ call max_all_cr(elemsize_max,elemsize_max_glob)
+
+ ! checks mesh
+ if( myrank == 0 ) then
+ if( distance_min_glob <= 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: GLL points minimum distance")
+ endif
+ if( distance_max_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: GLL points maximum distance")
+ endif
+ if( elemsize_min_glob <= 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: element minimum size")
+ endif
+ if( elemsize_max_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: element maximum size")
+ endif
+ endif
+
! min and max dimensions of the model
x_min = x_min_glob
x_max = x_max_glob
@@ -273,67 +331,25 @@
call min_all_cr(z_min,z_min_glob)
call max_all_cr(z_max,z_max_glob)
- ! element size
- elemsize_min = elemsize_min_glob
- elemsize_max = elemsize_max_glob
- call min_all_cr(elemsize_min,elemsize_min_glob)
- call max_all_cr(elemsize_max,elemsize_max_glob)
+ ! minimum period
+ pmax = pmax_glob
+ call max_all_cr(pmax,pmax_glob)
- ! checks mesh
- if( distance_min_glob <= 0.0_CUSTOM_REAL ) then
- call exit_mpi(myrank,"error: GLL points minimum distance")
+ ! time step
+ dt_suggested = dt_suggested_glob
+ call min_all_cr(dt_suggested,dt_suggested_glob)
+
+! determines global min/max values from all cpu partitions
+ if( DT_PRESENT ) then
+ ! courant number
+ cmax = cmax_glob
+ call max_all_cr(cmax,cmax_glob)
endif
- if( distance_max_glob >= HUGEVAL ) then
- call exit_mpi(myrank,"error: GLL points maximum distance")
- endif
- if( elemsize_min_glob <= 0.0_CUSTOM_REAL ) then
- call exit_mpi(myrank,"error: element minimum size")
- endif
- if( elemsize_max_glob >= HUGEVAL ) then
- call exit_mpi(myrank,"error: element maximum size")
- endif
-!! DK DK May 2009: added this to print the minimum and maximum number of elements
-!! DK DK May 2009: and points in the CUBIT + SCOTCH mesh
- call min_all_i(NSPEC_AB,NSPEC_AB_global_min)
- call max_all_i(NSPEC_AB,NSPEC_AB_global_max)
- call sum_all_i(NSPEC_AB,NSPEC_AB_global_sum)
+ !call world_size(sizeprocs)
- call min_all_i(NGLOB_AB,NGLOB_AB_global_min)
- call max_all_i(NGLOB_AB,NGLOB_AB_global_max)
- call sum_all_i(NGLOB_AB,NGLOB_AB_global_sum)
-
- call world_size(sizeprocs)
-
! outputs infos
if ( myrank == 0 ) then
- write(IMAIN,*)
- write(IMAIN,*) '********'
- write(IMAIN,*) 'minimum and maximum number of elements'
- write(IMAIN,*) 'and points in the CUBIT + SCOTCH mesh:'
- write(IMAIN,*)
- write(IMAIN,*) 'NSPEC_global_min = ',NSPEC_AB_global_min
- write(IMAIN,*) 'NSPEC_global_max = ',NSPEC_AB_global_max
- write(IMAIN,*) 'NSPEC_global_max / NSPEC_global_min imbalance = ',sngl(dble(NSPEC_AB_global_max) / dble(NSPEC_AB_global_min)),&
- ' = ',sngl((dble(NSPEC_AB_global_max) / dble(NSPEC_AB_global_min) - 1.d0) * 100.d0),' %'
- write(IMAIN,*) 'NSPEC_global_sum = ',NSPEC_AB_global_sum
- write(IMAIN,*)
- write(IMAIN,*) 'NGLOB_global_min = ',NGLOB_AB_global_min
- write(IMAIN,*) 'NGLOB_global_max = ',NGLOB_AB_global_max
- write(IMAIN,*) 'NGLOB_global_max / NGLOB_global_min imbalance = ',sngl(dble(NGLOB_AB_global_max) / dble(NGLOB_AB_global_min)),&
- ' = ',sngl((dble(NGLOB_AB_global_max) / dble(NGLOB_AB_global_min) - 1.d0) * 100.d0),' %'
- write(IMAIN,*) 'NGLOB_global_sum = ',NGLOB_AB_global_sum
- write(IMAIN,*)
- write(IMAIN,*) 'If you have elements of a single type (all acoustic, all elastic, all poroelastic, and without CPML)'
- write(IMAIN,*) 'in the whole mesh, then there should be no significant imbalance in the above numbers.'
- write(IMAIN,*) 'Otherwise, it is normal to have imbalance in elements and points because the domain decomposer'
- write(IMAIN,*) 'compensates for the different cost of different elements by partitioning them unevenly among processes.'
- write(IMAIN,*)
- write(IMAIN,*) '********'
- write(IMAIN,*) 'Model: P velocity min,max = ',vpmin_glob,vpmax_glob
- write(IMAIN,*) 'Model: S velocity min,max = ',vsmin_glob,vsmax_glob
- write(IMAIN,*) '********'
- write(IMAIN,*)
write(IMAIN,*) '*********************************************'
write(IMAIN,*) '*** Verification of simulation parameters ***'
write(IMAIN,*) '*********************************************'
@@ -358,6 +374,7 @@
write(IMAIN,*) '*** Max stability for wave velocities = ',cmax_glob
write(IMAIN,*)
endif
+ call flush_IMAIN()
endif
! returns the maximum velocity
@@ -381,6 +398,14 @@
! debug: for vtk output
if( SAVE_MESH_FILES ) then
call create_name_database(prname,myrank,LOCAL_PATH)
+
+ ! user output
+ if ( myrank == 0 ) then
+ write(IMAIN,*) 'saving VTK files for courant number and minimum period'
+ write(IMAIN,*)
+ call flush_IMAIN()
+ endif
+
! courant number
if( DT_PRESENT ) then
filename = trim(prname)//'res_courant_number'
@@ -393,6 +418,7 @@
call write_VTK_data_elem_cr(NSPEC_AB,NGLOB_AB, &
xstore,ystore,zstore,ibool, &
tmp2,filename)
+
deallocate(tmp1,tmp2)
endif
@@ -441,7 +467,7 @@
integer :: myrank
integer :: NSPEC_AB_global_min,NSPEC_AB_global_max,NSPEC_AB_global_sum
integer :: NGLOB_AB_global_min,NGLOB_AB_global_max,NGLOB_AB_global_sum
- integer :: ispec,sizeprocs
+ integer :: ispec !,sizeprocs
!********************************************************************************
@@ -620,24 +646,26 @@
call max_all_cr(vsmax,vsmax_glob)
! checks velocities
- if( vpmin_glob <= 0.0_CUSTOM_REAL ) then
- call exit_mpi(myrank,"error: vp minimum velocity")
+ if( myrank == 0 ) then
+ if( vpmin_glob <= 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: vp minimum velocity")
+ endif
+ if( vpmax_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: vp maximum velocity")
+ endif
+ if( vp2min_glob < 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: vp2 minimum velocity")
+ endif
+ if( vp2max_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: vp2 maximum velocity")
+ endif
+ if( vsmin_glob < 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: vs minimum velocity")
+ endif
+ if( vsmax_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: vs maximum velocity")
+ endif
endif
- if( vpmax_glob >= HUGEVAL ) then
- call exit_mpi(myrank,"error: vp maximum velocity")
- endif
- if( vp2min_glob < 0.0_CUSTOM_REAL ) then
- call exit_mpi(myrank,"error: vp2 minimum velocity")
- endif
- if( vp2max_glob >= HUGEVAL ) then
- call exit_mpi(myrank,"error: vp2 maximum velocity")
- endif
- if( vsmin_glob < 0.0_CUSTOM_REAL ) then
- call exit_mpi(myrank,"error: vs minimum velocity")
- endif
- if( vsmax_glob >= HUGEVAL ) then
- call exit_mpi(myrank,"error: vs maximum velocity")
- endif
! GLL point distance
distance_min = distance_min_glob
@@ -652,18 +680,20 @@
call max_all_cr(elemsize_max,elemsize_max_glob)
! checks mesh
- if( distance_min_glob <= 0.0_CUSTOM_REAL ) then
- call exit_mpi(myrank,"error: GLL points minimum distance")
+ if( myrank == 0 ) then
+ if( distance_min_glob <= 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: GLL points minimum distance")
+ endif
+ if( distance_max_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: GLL points maximum distance")
+ endif
+ if( elemsize_min_glob <= 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: element minimum size")
+ endif
+ if( elemsize_max_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: element maximum size")
+ endif
endif
- if( distance_max_glob >= HUGEVAL ) then
- call exit_mpi(myrank,"error: GLL points maximum distance")
- endif
- if( elemsize_min_glob <= 0.0_CUSTOM_REAL ) then
- call exit_mpi(myrank,"error: element minimum size")
- endif
- if( elemsize_max_glob >= HUGEVAL ) then
- call exit_mpi(myrank,"error: element maximum size")
- endif
!! DK DK May 2009: added this to print the minimum and maximum number of elements
!! DK DK May 2009: and points in the CUBIT + SCOTCH mesh
@@ -675,7 +705,7 @@
call max_all_i(NGLOB_AB,NGLOB_AB_global_max)
call sum_all_i(NGLOB_AB,NGLOB_AB_global_sum)
- call world_size(sizeprocs)
+ !call world_size(sizeprocs)
! outputs infos
if ( myrank == 0 ) then
@@ -717,6 +747,7 @@
write(IMAIN,*) '*** Max stability for wave velocities = ',cmax_glob
write(IMAIN,*)
endif
+ call flush_IMAIN()
endif
! returns the maximum velocity
@@ -956,9 +987,9 @@
do i=1,NGLLX-1
iglob_a = ibool(i ,j,k,ispec)
iglob_b = ibool(i+1,j,k,ispec)
- dist = sqrt( ( xstore(iglob_a) - xstore(iglob_b) )**2 &
+ dist = ( xstore(iglob_a) - xstore(iglob_b) )**2 &
+ ( ystore(iglob_a) - ystore(iglob_b) )**2 &
- + ( zstore(iglob_a) - zstore(iglob_b) )**2 )
+ + ( zstore(iglob_a) - zstore(iglob_b) )**2
if( dist < distance_min) distance_min = dist
if( dist > distance_max) distance_max = dist
enddo
@@ -971,9 +1002,9 @@
do j=1,NGLLY-1
iglob_a = ibool(i,j ,k,ispec)
iglob_b = ibool(i,j+1,k,ispec)
- dist = sqrt( ( xstore(iglob_a) - xstore(iglob_b) )**2 &
+ dist = ( xstore(iglob_a) - xstore(iglob_b) )**2 &
+ ( ystore(iglob_a) - ystore(iglob_b) )**2 &
- + ( zstore(iglob_a) - zstore(iglob_b) )**2 )
+ + ( zstore(iglob_a) - zstore(iglob_b) )**2
if( dist < distance_min) distance_min = dist
if( dist > distance_max) distance_max = dist
enddo
@@ -986,15 +1017,18 @@
do k=1,NGLLZ-1
iglob_a = ibool(i,j,k ,ispec)
iglob_b = ibool(i,j,k+1,ispec)
- dist = sqrt( ( xstore(iglob_a) - xstore(iglob_b) )**2 &
+ dist = ( xstore(iglob_a) - xstore(iglob_b) )**2 &
+ ( ystore(iglob_a) - ystore(iglob_b) )**2 &
- + ( zstore(iglob_a) - zstore(iglob_b) )**2 )
+ + ( zstore(iglob_a) - zstore(iglob_b) )**2
if( dist < distance_min) distance_min = dist
if( dist > distance_max) distance_max = dist
enddo
enddo
enddo
+ distance_min = sqrt( distance_min )
+ distance_max = sqrt( distance_max )
+
end subroutine get_GLL_minmaxdistance
!
@@ -1035,9 +1069,9 @@
do k = 1, NGLLZ, NGLLZ-1
ibool1 = ibool(i1,j,k,ispec)
ibool2 = ibool(i2,j,k,ispec)
- dist = sqrt((xstore(ibool1) - xstore(ibool2))**2 &
+ dist = (xstore(ibool1) - xstore(ibool2))**2 &
+ (ystore(ibool1) - ystore(ibool2))**2 &
- + (zstore(ibool1) - zstore(ibool2))**2)
+ + (zstore(ibool1) - zstore(ibool2))**2
if(dist < elemsize_min) elemsize_min = dist
if(dist > elemsize_max) elemsize_max = dist
enddo
@@ -1050,9 +1084,9 @@
do k = 1, NGLLZ, NGLLZ-1
ibool1 = ibool(i,j1,k,ispec)
ibool2 = ibool(i,j2,k,ispec)
- dist = sqrt((xstore(ibool1) - xstore(ibool2))**2 &
+ dist = (xstore(ibool1) - xstore(ibool2))**2 &
+ (ystore(ibool1) - ystore(ibool2))**2 &
- + (zstore(ibool1) - zstore(ibool2))**2)
+ + (zstore(ibool1) - zstore(ibool2))**2
if(dist < elemsize_min) elemsize_min = dist
if(dist > elemsize_max) elemsize_max = dist
enddo
@@ -1065,13 +1099,16 @@
do j = 1, NGLLY, NGLLY-1
ibool1 = ibool(i,j,k1,ispec)
ibool2 = ibool(i,j,k2,ispec)
- dist = sqrt((xstore(ibool1) - xstore(ibool2))**2 &
+ dist = (xstore(ibool1) - xstore(ibool2))**2 &
+ (ystore(ibool1) - ystore(ibool2))**2 &
- + (zstore(ibool1) - zstore(ibool2))**2)
+ + (zstore(ibool1) - zstore(ibool2))**2
if(dist < elemsize_min) elemsize_min = dist
if(dist > elemsize_max) elemsize_max = dist
enddo
enddo
+ elemsize_min = sqrt( elemsize_min )
+ elemsize_max = sqrt( elemsize_max )
+
end subroutine get_elem_minmaxsize
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90 2013-06-25 12:01:17 UTC (rev 22416)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90 2013-06-25 12:15:59 UTC (rev 22417)
@@ -340,6 +340,17 @@
stop 'error: please set USE_RICKER_TIME_FUNCTION to .true. in Par_file and recompile solver'
endif
+
+ ! absorbing conditions
+ ! stacey absorbing free surface must have also stacey turned on
+ if( STACEY_INSTEAD_OF_FREE_SURFACE ) STACEY_ABSORBING_CONDITIONS = .true.
+
+ ! only one absorbing condition is possible
+ if( PML_CONDITIONS ) then
+ if( STACEY_ABSORBING_CONDITIONS )&
+ stop 'error: please set STACEY_ABSORBING_CONDITIONS and STACEY_INSTEAD_OF_FREE_SURFACE to .false. in Par_file for PML'
+ endif
+
end subroutine read_parameter_file
!
More information about the CIG-COMMITS
mailing list