[cig-commits] r15688 - in seismo/3D/SPECFEM3D_SESAME/trunk: . EXAMPLES/homogeneous_halfspace EXAMPLES/layered_halfspace UTILS/ampuero_implicit_Clayton decompose_mesh_SCOTCH decompose_mesh_SCOTCH/OUTPUT_FILES
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Wed Sep 23 17:30:15 PDT 2009
Author: danielpeter
Date: 2009-09-23 17:30:14 -0700 (Wed, 23 Sep 2009)
New Revision: 15688
Modified:
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/Par_file
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/block_mesh.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/boundary_definition.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/cubit2specfem3d.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/boundary_definition.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/cubit2specfem3d.py
seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/ampuero_implicit_Clayton/ampuero_implicit_ABC_specfem3D.f90
seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90
seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/OUTPUT_FILES/boundary_definition.py
seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/OUTPUT_FILES/cubit2specfem3d.py
seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90
seismo/3D/SPECFEM3D_SESAME/trunk/initialize_simulation.f90
seismo/3D/SPECFEM3D_SESAME/trunk/iterate_time.f90
seismo/3D/SPECFEM3D_SESAME/trunk/prepare_timerun.f90
seismo/3D/SPECFEM3D_SESAME/trunk/read_mesh_databases.f90
seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90
seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D_par.f90
Log:
bug fixed in cubit2specfem3d.py; cleanup of code; new use of NSPEC_ATTENUATION_AB for attenuation array sizes
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/Par_file
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/Par_file 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/Par_file 2009-09-24 00:30:14 UTC (rev 15688)
@@ -38,7 +38,7 @@
SAVE_MESH_FILES = .true.
# path to store the local database file on each node
-LOCAL_PATH = /scratch/lustre/dpeter/SPECFEM3D_SESAME/DATABASES_MPI.block
+LOCAL_PATH = /scratch/lustre/dpeter/SPECFEM3D_SESAME/DATABASES_MPI.hom
# interval at which we output time step info and max of norm of displacement
NTSTEP_BETWEEN_OUTPUT_INFO = 500
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/block_mesh.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/block_mesh.py 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/block_mesh.py 2009-09-24 00:30:14 UTC (rev 15688)
@@ -13,7 +13,7 @@
# Meshing the volumes
-elementsize = 10000.0
+elementsize = 3750.0
cubit.cmd('volume 1 size '+str(elementsize))
cubit.cmd('mesh volume 1')
@@ -28,11 +28,11 @@
#### Define material properties for the 3 volumes ################
cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
cubit.cmd('block 1 attribute count 5')
-cubit.cmd('block 1 attribute index 1 1') # flag 1 for 1. material properties
-cubit.cmd('block 1 attribute index 2 2800')
-cubit.cmd('block 1 attribute index 3 1500')
-cubit.cmd('block 1 attribute index 4 2300')
-cubit.cmd('block 1 attribute index 5 100000.')
+cubit.cmd('block 1 attribute index 1 1') # flag for material properties: 1 for 1. volume
+cubit.cmd('block 1 attribute index 2 2800') # vp
+cubit.cmd('block 1 attribute index 3 1500') # vs
+cubit.cmd('block 1 attribute index 4 2300') # rho
+cubit.cmd('block 1 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
cubit.cmd('export mesh "top.e" dimension 3 overwrite')
cubit.cmd('save as "meshing.cub" overwrite')
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/boundary_definition.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/boundary_definition.py 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/boundary_definition.py 2009-09-24 00:30:14 UTC (rev 15688)
@@ -86,7 +86,7 @@
# if dz <= 0.001 and dn < 0.2:
# top_surf.append(k)
-#daniel
+ #box lengths
x_len = abs( xmax_box - xmin_box)
y_len = abs( ymax_box - ymin_box)
z_len = abs( zmax_box - zmin_box)
@@ -280,8 +280,6 @@
id_nodeset=cubit.get_next_nodeset_id()
id_block=cubit.get_next_block_id()
- #daniel
- #print "##next nodeset: "+ str(id_nodeset) + " block: " + str(id_block)
if obj == 'hex':
txt='hex in node in surface'
@@ -300,7 +298,7 @@
txt2 = "block "+str(id_block)+" name '"+name+"'"
else:
txt1=''
- #daniel: do not execute as block id might be wrong
+ # do not execute: block id might be wrong
print "##block "+str(id_block)+" name '"+name+"_notsupported (only hex,face,edge,node)'"
txt2=''
@@ -312,7 +310,6 @@
parallel=keys.get('parallel',True)
closed=keys.get('closed',False)
if not closed:
- #dp
print "##open region"
if parallel:
@@ -324,7 +321,6 @@
entities=args[0]
print entities
for entity in entities:
- #dp
print "##entity: "+str(entity)
build_block_side(topo,entity+'_topo',obj=entity)
@@ -336,7 +332,6 @@
build_block_side(ymax,entity+'_abs_ymax',obj=entity)
build_block_side(bottom,entity+'_abs_bottom',obj=entity)
else:
- #dp
print "##closed region"
surf=define_absorbing_surf_sphere()
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/cubit2specfem3d.py 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/cubit2specfem3d.py 2009-09-24 00:30:14 UTC (rev 15688)
@@ -371,9 +371,17 @@
if nattrib >= 3:
vs=cubit.get_block_attribute_value(block,2)
if nattrib >= 4:
- rho=vs=cubit.get_block_attribute_value(block,3)
+ #density
+ rho=cubit.get_block_attribute_value(block,3)
if nattrib == 5:
- q=vs=cubit.get_block_attribute_value(block,4)
+ #Q_flag
+ q=cubit.get_block_attribute_value(block,4)
+ # for q to be valid: it must be an integer flag between 1 and 13
+ # (see constants.h for IATTENUATION_SEDIMENT_40, etc. )
+ if q < 0 or q > 13:
+ print 'error, q flag invalid:', q
+ print ' check with constants.h for IATTENUATION flags'
+ break
elif flag < 0:
vel=name
attrib=cubit.get_block_attribute_value(block,1)
@@ -447,7 +455,8 @@
rho=(1.6612*vp-0.472*vp**2+0.0671*vp**3-0.0043*vp**4+0.000106*vp**4)*m2km
txt='%3i %20f %20f %20f %1i %1i\n' % (properties[0],rho,vel,vel/(3**.5),0,0)
elif type(vel) != str:
- txt='%3i %20f %20f %20f %1i %1i\n' % (properties[0],properties[3],properties[1],properties[2],0,0)
+ #format nummaterials file: #material_id #rho #vp #vs #Q_flag #0
+ txt='%3i %20f %20f %20f %2i %1i\n' % (properties[0],properties[3],properties[1],properties[2],properties[4],0)
else:
txt='%3i %s \n' % (properties[0],properties[1])
elif flag < 0:
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py 2009-09-24 00:30:14 UTC (rev 15688)
@@ -71,25 +71,25 @@
#### Define material properties for the 3 volumes ################
cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
cubit.cmd('block 1 attribute count 5')
-cubit.cmd('block 1 attribute index 1 1 ')
-cubit.cmd('block 1 attribute index 2 2800 ')
-cubit.cmd('block 1 attribute index 3 1500 ')
-cubit.cmd('block 1 attribute index 4 2300 ')
-cubit.cmd('block 1 attribute index 5 100000. ')
+cubit.cmd('block 1 attribute index 1 1 ') # volume 1
+cubit.cmd('block 1 attribute index 2 2800 ') # vp
+cubit.cmd('block 1 attribute index 3 1500 ') # vs
+cubit.cmd('block 1 attribute index 4 2300 ') # rho
+cubit.cmd('block 1 attribute index 5 1 ') # Q_flag
cubit.cmd('block 2 attribute count 5')
-cubit.cmd('block 2 attribute index 1 2 ')
+cubit.cmd('block 2 attribute index 1 2 ') # volume 2
cubit.cmd('block 2 attribute index 2 7500 ')
cubit.cmd('block 2 attribute index 3 4300 ')
cubit.cmd('block 2 attribute index 4 3200 ')
-cubit.cmd('block 2 attribute index 5 100000. ')
+cubit.cmd('block 2 attribute index 5 1')
cubit.cmd('block 3 attribute count 5')
-cubit.cmd('block 3 attribute index 1 2 ')
+cubit.cmd('block 3 attribute index 1 2 ') # same material properties as for volume 2
cubit.cmd('block 3 attribute index 2 7500 ')
cubit.cmd('block 3 attribute index 3 4300 ')
cubit.cmd('block 3 attribute index 4 3200 ')
-cubit.cmd('block 3 attribute index 5 100000. ')
+cubit.cmd('block 3 attribute index 5 1')
cubit.cmd('export mesh "top.e" dimension 3 overwrite')
cubit.cmd('save as "meshing.cub" overwrite')
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py 2009-09-24 00:30:14 UTC (rev 15688)
@@ -71,25 +71,25 @@
#### Define material properties for the 3 volumes ################
cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
cubit.cmd('block 1 attribute count 5')
-cubit.cmd('block 1 attribute index 1 1 ')
-cubit.cmd('block 1 attribute index 2 2800 ')
-cubit.cmd('block 1 attribute index 3 1500 ')
-cubit.cmd('block 1 attribute index 4 2300 ')
-cubit.cmd('block 1 attribute index 5 100000. ')
+cubit.cmd('block 1 attribute index 1 1 ') # volume 1
+cubit.cmd('block 1 attribute index 2 2800 ') # vp
+cubit.cmd('block 1 attribute index 3 1500 ') # vs
+cubit.cmd('block 1 attribute index 4 2300 ') # rho
+cubit.cmd('block 1 attribute index 5 1 ') # Q_flag
cubit.cmd('block 2 attribute count 5')
-cubit.cmd('block 2 attribute index 1 2 ')
+cubit.cmd('block 2 attribute index 1 2 ') # volume 2
cubit.cmd('block 2 attribute index 2 7500 ')
cubit.cmd('block 2 attribute index 3 4300 ')
cubit.cmd('block 2 attribute index 4 3200 ')
-cubit.cmd('block 2 attribute index 5 100000. ')
+cubit.cmd('block 2 attribute index 5 1 ')
cubit.cmd('block 3 attribute count 5')
-cubit.cmd('block 3 attribute index 1 2 ')
+cubit.cmd('block 3 attribute index 1 2 ') # same properties as for volume 2
cubit.cmd('block 3 attribute index 2 7500 ')
cubit.cmd('block 3 attribute index 3 4300 ')
cubit.cmd('block 3 attribute index 4 3200 ')
-cubit.cmd('block 3 attribute index 5 100000. ')
+cubit.cmd('block 3 attribute index 5 1 ')
cubit.cmd('export mesh "top.e" dimension 3 overwrite')
cubit.cmd('save as "meshing.cub" overwrite')
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/boundary_definition.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/boundary_definition.py 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/boundary_definition.py 2009-09-24 00:30:14 UTC (rev 15688)
@@ -86,7 +86,7 @@
# if dz <= 0.001 and dn < 0.2:
# top_surf.append(k)
-#daniel
+ #box lengths
x_len = abs( xmax_box - xmin_box)
y_len = abs( ymax_box - ymin_box)
z_len = abs( zmax_box - zmin_box)
@@ -280,8 +280,6 @@
id_nodeset=cubit.get_next_nodeset_id()
id_block=cubit.get_next_block_id()
- #daniel
- #print "##next nodeset: "+ str(id_nodeset) + " block: " + str(id_block)
if obj == 'hex':
txt='hex in node in surface'
@@ -300,7 +298,7 @@
txt2 = "block "+str(id_block)+" name '"+name+"'"
else:
txt1=''
- #daniel: do not execute as block id might be wrong
+ # do not execute: block id might be wrong
print "##block "+str(id_block)+" name '"+name+"_notsupported (only hex,face,edge,node)'"
txt2=''
@@ -312,7 +310,6 @@
parallel=keys.get('parallel',True)
closed=keys.get('closed',False)
if not closed:
- #dp
print "##open region"
if parallel:
@@ -324,7 +321,6 @@
entities=args[0]
print entities
for entity in entities:
- #dp
print "##entity: "+str(entity)
build_block_side(topo,entity+'_topo',obj=entity)
@@ -336,7 +332,6 @@
build_block_side(ymax,entity+'_abs_ymax',obj=entity)
build_block_side(bottom,entity+'_abs_bottom',obj=entity)
else:
- #dp
print "##closed region"
surf=define_absorbing_surf_sphere()
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/cubit2specfem3d.py 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/cubit2specfem3d.py 2009-09-24 00:30:14 UTC (rev 15688)
@@ -371,9 +371,17 @@
if nattrib >= 3:
vs=cubit.get_block_attribute_value(block,2)
if nattrib >= 4:
- rho=vs=cubit.get_block_attribute_value(block,3)
+ #density
+ rho=cubit.get_block_attribute_value(block,3)
if nattrib == 5:
- q=vs=cubit.get_block_attribute_value(block,4)
+ #Q_flag
+ q=cubit.get_block_attribute_value(block,4)
+ # for q to be valid: it must be an integer flag between 1 and 13
+ # (see constants.h for IATTENUATION_SEDIMENT_40, etc. )
+ if q < 0 or q > 13:
+ print 'error, q flag invalid:', q
+ print ' check with constants.h for IATTENUATION flags'
+ break
elif flag < 0:
vel=name
attrib=cubit.get_block_attribute_value(block,1)
@@ -447,7 +455,8 @@
rho=(1.6612*vp-0.472*vp**2+0.0671*vp**3-0.0043*vp**4+0.000106*vp**4)*m2km
txt='%3i %20f %20f %20f %1i %1i\n' % (properties[0],rho,vel,vel/(3**.5),0,0)
elif type(vel) != str:
- txt='%3i %20f %20f %20f %1i %1i\n' % (properties[0],properties[3],properties[1],properties[2],0,0)
+ #format nummaterials file: #material_id #rho #vp #vs #Q_flag #0
+ txt='%3i %20f %20f %20f %2i %1i\n' % (properties[0],properties[3],properties[1],properties[2],properties[4],0)
else:
txt='%3i %s \n' % (properties[0],properties[1])
elif flag < 0:
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in 2009-09-24 00:30:14 UTC (rev 15688)
@@ -114,7 +114,6 @@
# solver objects with statically allocated arrays; dependent upon
# values_from_mesher.h
-#daniel: added files
SOLVER_ARRAY_OBJECTS = \
$O/assemble_MPI_scalar.o \
$O/assemble_MPI_vector.o \
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/ampuero_implicit_Clayton/ampuero_implicit_ABC_specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/ampuero_implicit_Clayton/ampuero_implicit_ABC_specfem3D.f90 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/ampuero_implicit_Clayton/ampuero_implicit_ABC_specfem3D.f90 2009-09-24 00:30:14 UTC (rev 15688)
@@ -1186,7 +1186,7 @@
duzdyl_plus_duydzl = duzdyl + duydzl
! precompute terms for attenuation if needed
- if(ATTENUATION_VAL) then
+ if(ATTENUATION) then
! compute deviatoric strain
epsilon_trace_over_3 = ONE_THIRD * (duxdxl + duydyl + duzdzl)
@@ -1207,7 +1207,7 @@
mul = mustore(i,j,k,ispec)
! use unrelaxed parameters if attenuation
- if(ATTENUATION_VAL) mul = mul * one_minus_sum_beta_use
+ if(ATTENUATION) mul = mul * one_minus_sum_beta_use
lambdalplus2mul = kappal + FOUR_THIRDS * mul
lambdal = lambdalplus2mul - 2.*mul
@@ -1222,7 +1222,7 @@
sigma_yz = mul*duzdyl_plus_duydzl
! subtract memory variables if attenuation
- if(ATTENUATION_VAL) then
+ if(ATTENUATION) then
do i_sls = 1,N_SLS
R_xx_val = R_xx(i,j,k,ispec,i_sls)
R_yy_val = R_yy(i,j,k,ispec,i_sls)
@@ -1303,7 +1303,7 @@
! update memory variables based upon the Runge-Kutta scheme
- if(ATTENUATION_VAL) then
+ if(ATTENUATION) then
! use Runge-Kutta scheme to march in time
do i_sls = 1,N_SLS
@@ -1351,7 +1351,7 @@
enddo
! save deviatoric strain for Runge-Kutta scheme
- if(ATTENUATION_VAL) then
+ if(ATTENUATION) then
epsilondev_xx(:,:,:,ispec) = epsilondev_xx_loc(:,:,:)
epsilondev_yy(:,:,:,ispec) = epsilondev_yy_loc(:,:,:)
epsilondev_xy(:,:,:,ispec) = epsilondev_xy_loc(:,:,:)
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90 2009-09-24 00:30:14 UTC (rev 15688)
@@ -23,13 +23,13 @@
!
!=====================================================================
-subroutine compute_forces_with_Deville(NSPEC_AB,NGLOB_AB,ATTENUATION_VAL,displ,accel, &
+subroutine compute_forces_with_Deville(NSPEC_AB,NGLOB_AB,ATTENUATION,displ,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source, &
hdur,hdur_gaussian,t_cmt,dt,stf,t0,sourcearrays, & !pll
- one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,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,iflag_attenuation_store, &
ABSORBING_CONDITIONS, &
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2DMAX_XMIN_XMAX_ext,NSPEC2DMAX_YMIN_YMAX_ext,&
@@ -142,13 +142,14 @@
epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
real(kind=CUSTOM_REAL) epsilon_trace_over_3
- logical :: ATTENUATION_VAL
+ logical :: ATTENUATION
+ integer :: NSPEC_ATTENUATION_AB
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: iflag_attenuation_store
real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION) :: one_minus_sum_beta
real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: factor_common, alphaval,betaval,gammaval
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB,N_SLS) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS) :: &
R_xx,R_yy,R_xy,R_xz,R_yz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB) :: &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
! Stacey conditions
@@ -309,7 +310,7 @@
kappal = kappastore(i,j,k,ispec)
mul = mustore(i,j,k,ispec)
- if(ATTENUATION_VAL) then
+ if(ATTENUATION) then
! compute deviatoric strain
epsilon_trace_over_3 = ONE_THIRD * (duxdxl + duydyl + duzdzl)
epsilondev_xx_loc(i,j,k) = duxdxl - epsilon_trace_over_3
@@ -335,7 +336,7 @@
sigma_yz = mul*duzdyl_plus_duydzl
! subtract memory variables if attenuation
- if(ATTENUATION_VAL) then
+ if(ATTENUATION) then
do i_sls = 1,N_SLS
R_xx_val = R_xx(i,j,k,ispec,i_sls)
R_yy_val = R_yy(i,j,k,ispec,i_sls)
@@ -347,7 +348,9 @@
sigma_yz = sigma_yz - R_yz(i,j,k,ispec,i_sls)
enddo
endif
+
+
! form dot product with test vector, symmetric form
tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
@@ -456,7 +459,7 @@
accel(3,iglob) = accel(3,iglob) - fac1*newtempz1(i,j,k) - fac2*newtempz2(i,j,k) - fac3*newtempz3(i,j,k)
! update memory variables based upon the Runge-Kutta scheme
- if(ATTENUATION_VAL) then
+ if(ATTENUATION) then
! use Runge-Kutta scheme to march in time
do i_sls = 1,N_SLS
@@ -504,7 +507,7 @@
enddo
! save deviatoric strain for Runge-Kutta scheme
- if(ATTENUATION_VAL) then
+ if(ATTENUATION) then
epsilondev_xx(:,:,:,ispec) = epsilondev_xx_loc(:,:,:)
epsilondev_yy(:,:,:,ispec) = epsilondev_yy_loc(:,:,:)
epsilondev_xy(:,:,:,ispec) = epsilondev_xy_loc(:,:,:)
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90 2009-09-24 00:30:14 UTC (rev 15688)
@@ -339,6 +339,7 @@
do i = 1, NGLLX
! check if the material is known or unknown
if (mat_ext_mesh(1,ispec) > 0) then
+ ! materials_ext_mesh format: #index1 = rho #index2 = vp #index3 = vs #index4 = Q_flag #index5 = 0
rhostore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(1,ispec))
vpstore(i,j,k,ispec) = materials_ext_mesh(2,mat_ext_mesh(1,ispec))
vsstore(i,j,k,ispec) = materials_ext_mesh(3,mat_ext_mesh(1,ispec))
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/OUTPUT_FILES/boundary_definition.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/OUTPUT_FILES/boundary_definition.py 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/OUTPUT_FILES/boundary_definition.py 2009-09-24 00:30:14 UTC (rev 15688)
@@ -49,6 +49,8 @@
absorbing_surf_ymax=[]
absorbing_surf_bottom=[]
top_surf=[]
+
+
list_vol=cubit.parse_cubit_list("volume","all")
init_n_vol=len(list_vol)
zmax_box=cubit.get_total_bounding_box("volume",list_vol)[7]
@@ -58,31 +60,73 @@
ymin_box=cubit.get_total_bounding_box("volume",list_vol)[3]
ymax_box=cubit.get_total_bounding_box("volume",list_vol)[4]
list_surf=cubit.parse_cubit_list("surface","all")
+# for k in list_surf:
+# center_point = cubit.get_center_point("surface", k)
+# if abs((center_point[0] - xmin_box)/xmin_box) <= 0.005:
+# absorbing_surf_xmin.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[0] - xmax_box)/xmax_box) <= 0.005:
+# absorbing_surf_xmax.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[1] - ymin_box)/ymin_box) <= 0.005:
+# absorbing_surf_ymin.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[1] - ymax_box)/ymax_box) <= 0.005:
+# absorbing_surf_ymax.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[2] - zmin_box)/zmin_box) <= 0.005:
+# absorbing_surf_bottom.append(k)
+# absorbing_surf.append(k)
+# else:
+# sbox=cubit.get_bounding_box('surface',k)
+# dz=abs((sbox[7] - zmax_box)/zmax_box)
+# normal=cubit.get_surface_normal(k)
+# zn=normal[2]
+# dn=abs(zn-1)
+# if dz <= 0.001 and dn < 0.2:
+# top_surf.append(k)
+
+ #box lengths
+ x_len = abs( xmax_box - xmin_box)
+ y_len = abs( ymax_box - ymin_box)
+ z_len = abs( zmax_box - zmin_box)
+
+ print '##boundary box: '
+ print '## x length: ' + str(x_len)
+ print '## y length: ' + str(y_len)
+ print '## z length: ' + str(z_len)
+
+ # tolerance parameters
+ absorbing_surface_distance_tolerance=0.005
+ topographic_surface_distance_tolerance=0.001
+ topographic_surface_normal_tolerance=0.2
+
for k in list_surf:
center_point = cubit.get_center_point("surface", k)
- if abs((center_point[0] - xmin_box)/xmin_box) <= 0.005:
+ if abs((center_point[0] - xmin_box)/x_len) <= absorbing_surface_distance_tolerance:
absorbing_surf_xmin.append(k)
absorbing_surf.append(k)
- elif abs((center_point[0] - xmax_box)/xmax_box) <= 0.005:
+ elif abs((center_point[0] - xmax_box)/x_len) <= absorbing_surface_distance_tolerance:
absorbing_surf_xmax.append(k)
absorbing_surf.append(k)
- elif abs((center_point[1] - ymin_box)/ymin_box) <= 0.005:
+ elif abs((center_point[1] - ymin_box)/y_len) <= absorbing_surface_distance_tolerance:
absorbing_surf_ymin.append(k)
absorbing_surf.append(k)
- elif abs((center_point[1] - ymax_box)/ymax_box) <= 0.005:
+ elif abs((center_point[1] - ymax_box)/y_len) <= absorbing_surface_distance_tolerance:
absorbing_surf_ymax.append(k)
absorbing_surf.append(k)
- elif abs((center_point[2] - zmin_box)/zmin_box) <= 0.005:
+ elif abs((center_point[2] - zmin_box)/z_len) <= absorbing_surface_distance_tolerance:
absorbing_surf_bottom.append(k)
absorbing_surf.append(k)
else:
sbox=cubit.get_bounding_box('surface',k)
- dz=abs((sbox[7] - zmax_box)/zmax_box)
+ dz=abs((sbox[7] - zmax_box)/z_len)
normal=cubit.get_surface_normal(k)
zn=normal[2]
dn=abs(zn-1)
- if dz <= 0.001 and dn < 0.2:
+ if dz <= topographic_surface_distance_tolerance and dn < topographic_surface_normal_tolerance:
top_surf.append(k)
+
return absorbing_surf,absorbing_surf_xmin,absorbing_surf_xmax,absorbing_surf_ymin,absorbing_surf_ymax,absorbing_surf_bottom,top_surf
def define_absorbing_surf_nopar():
@@ -235,6 +279,8 @@
sys.exit()
id_nodeset=cubit.get_next_nodeset_id()
id_block=cubit.get_next_block_id()
+
+
if obj == 'hex':
txt='hex in node in surface'
txt1='block '+str(id_block)+ ' '+ txt +' '+str(list(surf_list))
@@ -252,7 +298,11 @@
txt2 = "block "+str(id_block)+" name '"+name+"'"
else:
txt1=''
- txt2="block "+str(id_block)+" name '"+name+"_notsupported (only hex,face,edge,node)'"
+ # do not execute: block id might be wrong
+ print "##block "+str(id_block)+" name '"+name+"_notsupported (only hex,face,edge,node)'"
+ txt2=''
+
+
cubit.cmd(txt1)
cubit.cmd(txt2)
@@ -260,6 +310,8 @@
parallel=keys.get('parallel',True)
closed=keys.get('closed',False)
if not closed:
+ print "##open region"
+
if parallel:
surf,xmin,xmax,ymin,ymax,bottom,topo=define_absorbing_surf()
else:
@@ -269,6 +321,8 @@
entities=args[0]
print entities
for entity in entities:
+ print "##entity: "+str(entity)
+
build_block_side(topo,entity+'_topo',obj=entity)
build_block_side(surf,entity+'_abs',obj=entity)
if parallel:
@@ -278,6 +332,8 @@
build_block_side(ymax,entity+'_abs_ymax',obj=entity)
build_block_side(bottom,entity+'_abs_bottom',obj=entity)
else:
+ print "##closed region"
+
surf=define_absorbing_surf_sphere()
v_list,name_list=define_block()
build_block(v_list,name_list)
@@ -291,9 +347,12 @@
#define_bc(entities,parallel=True)
#define_bc(entities,parallel=False)
#define_bc(entities,parallel=False,closed=True)
-entities=['surface','face']
-define_bc(entities,parallel=True)
+
+# call
+#entities=['surface','face']
+#define_bc(entities,parallel=True)
+
#block 1 attribute count 5
#block 2 attribute count 0
#block 2 name '/prova/interface1'
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/OUTPUT_FILES/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/OUTPUT_FILES/cubit2specfem3d.py 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/OUTPUT_FILES/cubit2specfem3d.py 2009-09-24 00:30:14 UTC (rev 15688)
@@ -104,9 +104,9 @@
#############################################################################
+import cubit
-
class mtools(object):
"""docstring for ciao"""
def __init__(self,frequency,list_surf,list_vp):
@@ -371,9 +371,17 @@
if nattrib >= 3:
vs=cubit.get_block_attribute_value(block,2)
if nattrib >= 4:
- rho=vs=cubit.get_block_attribute_value(block,3)
+ #density
+ rho=cubit.get_block_attribute_value(block,3)
if nattrib == 5:
- q=vs=cubit.get_block_attribute_value(block,4)
+ #Q_flag
+ q=cubit.get_block_attribute_value(block,4)
+ # for q to be valid: it must be an integer flag between 1 and 13
+ # (see constants.h for IATTENUATION_SEDIMENT_40, etc. )
+ if q < 0 or q > 13:
+ print 'error, q flag invalid:', q
+ print ' check with constants.h for IATTENUATION flags'
+ break
elif flag < 0:
vel=name
attrib=cubit.get_block_attribute_value(block,1)
@@ -447,7 +455,8 @@
rho=(1.6612*vp-0.472*vp**2+0.0671*vp**3-0.0043*vp**4+0.000106*vp**4)*m2km
txt='%3i %20f %20f %20f %1i %1i\n' % (properties[0],rho,vel,vel/(3**.5),0,0)
elif type(vel) != str:
- txt='%3i %20f %20f %20f %1i %1i\n' % (properties[0],properties[3],properties[1],properties[2],0,0)
+ #format nummaterials file: #material_id #rho #vp #vs #Q_flag #0
+ txt='%3i %20f %20f %20f %2i %1i\n' % (properties[0],properties[3],properties[1],properties[2],properties[4],0)
else:
txt='%3i %s \n' % (properties[0],properties[1])
elif flag < 0:
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2009-09-24 00:30:14 UTC (rev 15688)
@@ -78,6 +78,7 @@
if((num_elmnt > nspec) .or. (num_elmnt < 1) ) stop "ERROR : Invalid mesh file."
end do
close(98)
+ print*, 'total number of spectral elements:'
print*, 'nspec = ', nspec
open(unit=98, file='./OUTPUT_FILES/materials_file', status='old', form='formatted')
@@ -98,12 +99,14 @@
!if(num_node /= inode) stop "ERROR : Invalid nodes_coords file."
end do
close(98)
+ print*, 'total number of nodes: '
print*, 'nnodes = ', nnodes
count_def_mat = 0
count_undef_mat = 0
open(unit=98, file='./OUTPUT_FILES/nummaterial_velocity_file', status='old', form='formatted')
read(98,*,iostat=ierr) num_mat
+ print *,'materials:'
do while (ierr == 0)
print*, 'num_mat = ',num_mat
if(num_mat /= -1) then
@@ -115,14 +118,21 @@
end do
close(98)
- print*, count_def_mat, count_undef_mat
+ print*, 'defined = ',count_def_mat, 'undefined = ',count_undef_mat
allocate(mat_prop(5,count_def_mat))
allocate(undef_mat_prop(5,count_undef_mat))
open(unit=98, file='./OUTPUT_FILES/nummaterial_velocity_file', status='old', form='formatted')
do imat=1,count_def_mat
+ ! format:# material_id # rho # vp # vs # Q_flag # 0
read(98,*) num_mat, mat_prop(1,num_mat),mat_prop(2,num_mat),mat_prop(3,num_mat),mat_prop(4,num_mat),mat_prop(5,num_mat)
if(num_mat < 0 .or. num_mat > count_def_mat) stop "ERROR : Invalid nummaterial_velocity_file file."
+
+ !checks attenuation flag with integer range as defined in constants.h like IATTENUATION_SEDIMENTS_40, ....
+ if( int(mat_prop(4,num_mat)) > 13 ) then
+ stop 'wrong attenuation flag in mesh: too large, not supported yet - check with constants.h'
+ endif
+
end do
do imat=1,count_undef_mat
read(98,'(5A30)') undef_mat_prop(1,imat),undef_mat_prop(2,imat),undef_mat_prop(3,imat),undef_mat_prop(4,imat), &
@@ -140,6 +150,7 @@
nodes_ibelm_xmin(3,ispec2D), nodes_ibelm_xmin(4,ispec2D)
end do
close(98)
+ print*, 'absorbing boundaries:'
print*, 'nspec2D_xmin = ', nspec2D_xmin
open(unit=98, file='./OUTPUT_FILES/absorbing_surface_file_xmax', status='old', form='formatted')
@@ -207,7 +218,8 @@
used_nodes_elmnts(elmnts(inode,ispec)) = used_nodes_elmnts(elmnts(inode,ispec)) + 1
enddo
enddo
- print *, minval(used_nodes_elmnts(:)), maxval(used_nodes_elmnts(:))
+ print *, 'nodes valence: '
+ print *, 'min = ',minval(used_nodes_elmnts(:)),'max = ', maxval(used_nodes_elmnts(:))
do inode = 1, nnodes
if (.not. mask_nodes_elmnts(inode)) then
stop 'ERROR : nodes not used.'
@@ -226,6 +238,7 @@
call mesh2dual_ncommonnodes(nspec, nnodes, nsize, sup_neighbour, elmnts, xadj, adjncy, nnodes_elmnts, &
nodes_elmnts, max_neighbour, 1)
+ print*, 'mesh2dual: '
print*, 'max_neighbour = ',max_neighbour
nb_edges = xadj(nspec+1)
@@ -328,7 +341,11 @@
close(15)
end do
-
+ print*, 'partitions: '
+ print*, 'num = ',nparts
+ print*, 'files in directory: OUTPUT_FILES/'
+ print*
+ print*, 'finished successfully'
end program pre_meshfem3D
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2009-09-24 00:30:14 UTC (rev 15688)
@@ -421,6 +421,7 @@
write(IIN_database,*) count_def_mat,count_undef_mat
do i = 1, count_def_mat
+ ! format: # rho # vp # vs # Q_flag # 0
write(IIN_database,*) mat_prop(1,i), mat_prop(2,i), mat_prop(3,i), mat_prop(4,i), mat_prop(5,i)
end do
do i = 1, count_undef_mat
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90 2009-09-24 00:30:14 UTC (rev 15688)
@@ -474,6 +474,7 @@
allocate(materials_ext_mesh(5,nmat_ext_mesh))
allocate(undef_mat_prop(5,nundefMat_ext_mesh))
do imat = 1, nmat_ext_mesh
+ ! format: # rho # vp # vs # Q_flag # 0
read(IIN,*) materials_ext_mesh(1,imat), materials_ext_mesh(2,imat), materials_ext_mesh(3,imat), &
materials_ext_mesh(4,imat), materials_ext_mesh(5,imat)
end do
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/initialize_simulation.f90 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/initialize_simulation.f90 2009-09-24 00:30:14 UTC (rev 15688)
@@ -79,10 +79,16 @@
read(27) NSPEC_AB
read(27) NGLOB_AB
- !pll
- NSPEC_ATTENUATION_AB = NSPEC_AB
close(27)
+ if( ATTENUATION ) then
+ !pll
+ NSPEC_ATTENUATION_AB = NSPEC_AB
+ else
+ ! if attenuation is off, set dummy size of arrays to one
+ NSPEC_ATTENUATION_AB = 1
+ endif
+
! open main output file, only written to by process 0
if(myrank == 0 .and. IMAIN /= ISTANDARD_OUTPUT) &
open(unit=IMAIN,file=trim(OUTPUT_FILES)//'/output_solver.txt',status='unknown')
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/iterate_time.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/iterate_time.f90 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/iterate_time.f90 2009-09-24 00:30:14 UTC (rev 15688)
@@ -58,21 +58,7 @@
! *********************************************************
do it = 1,NSTEP
-
-
-!check stability
- do i=1,3
- Usolidnorm = maxval(abs(displ(i,:)))
- Usolidnorm_index = maxloc(abs(displ(i,:)))
- if(Usolidnorm > 1.e+15 ) then
- print*,' stability issue:',myrank
- print*,' norm: ',Usolidnorm,displ(i,Usolidnorm_index(1)),i
- print*,' index: ',Usolidnorm_index(1)
- print*,' x/y/z: ',xstore(Usolidnorm_index(1)),ystore(Usolidnorm_index(1)),zstore(Usolidnorm_index(1))
- print*,' time step: ',it
- call exit_MPI(myrank,'forward simulation became unstable and blew up')
- endif
- enddo
+
! compute the maximum of the norm of the displacement
! in all the slices using an MPI reduction
! and output timestamp file to check that simulation is running fine
@@ -107,8 +93,40 @@
write(IMAIN,*) 'Max norm displacement vector U in all slices (m) = ',Usolidnorm_all
! if (SIMULATION_TYPE == 3) write(IMAIN,*) &
! 'Max norm displacement vector U (backward) in all slices (m) = ',b_Usolidnorm_all
+
+! compute estimated remaining simulation time
+ t_remain = (NSTEP - it) * (tCPU/dble(it))
+ int_t_remain = int(t_remain)
+ ihours_remain = int_t_remain / 3600
+ iminutes_remain = (int_t_remain - 3600*ihours_remain) / 60
+ iseconds_remain = int_t_remain - 3600*ihours_remain - 60*iminutes_remain
+ write(IMAIN,*) 'Time steps done = ',it,' out of ',NSTEP
+ write(IMAIN,*) 'Time steps remaining = ',NSTEP - it
+ write(IMAIN,*) 'Estimated remaining time in seconds = ',t_remain
+ write(IMAIN,"(' Estimated remaining time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") &
+ ihours_remain,iminutes_remain,iseconds_remain
+
+! compute estimated total simulation time
+ t_total = t_remain + tCPU
+ int_t_total = int(t_total)
+ ihours_total = int_t_total / 3600
+ iminutes_total = (int_t_total - 3600*ihours_total) / 60
+ iseconds_total = int_t_total - 3600*ihours_total - 60*iminutes_total
+ write(IMAIN,*) 'Estimated total run time in seconds = ',t_total
+ write(IMAIN,"(' Estimated total run time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") &
+ ihours_total,iminutes_total,iseconds_total
+ write(IMAIN,*) 'We have done ',sngl(100.d0*dble(it)/dble(NSTEP)),'% of that'
+
+ if(it < 100) then
+ write(IMAIN,*) '************************************************************'
+ write(IMAIN,*) '**** BEWARE: the above time estimates are not reliable'
+ write(IMAIN,*) '**** because fewer than 100 iterations have been performed'
+ write(IMAIN,*) '************************************************************'
+ endif
write(IMAIN,*)
+
+
! write time stamp file to give information about progression of simulation
write(outputname,"('/timestamp',i6.6)") it
open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='unknown')
@@ -163,7 +181,7 @@
kappastore,mustore,jacobian,ibool,ispec_is_inner_ext_mesh,.false., &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source, &
hdur,hdur_gaussian,t_cmt,dt,stf,t0,sourcearrays, &
- one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,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,iflag_attenuation_store,ABSORBING_CONDITIONS, &
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2DMAX_XMIN_XMAX_ext,NSPEC2DMAX_YMIN_YMAX_ext, &
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom, &
@@ -189,7 +207,7 @@
kappastore,mustore,jacobian,ibool,ispec_is_inner_ext_mesh,.true., &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source, &
hdur,hdur_gaussian,t_cmt,dt,stf,t0,sourcearrays, &
- one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,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,iflag_attenuation_store,ABSORBING_CONDITIONS, &
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2DMAX_XMIN_XMAX_ext,NSPEC2DMAX_YMIN_YMAX_ext, &
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom, &
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/prepare_timerun.f90 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/prepare_timerun.f90 2009-09-24 00:30:14 UTC (rev 15688)
@@ -298,6 +298,8 @@
!pll, to put elsewhere
+ ! note: currently, they need to be defined here, as they are used in the routine arguments
+ ! for compute_forces_with_Deville()
allocate(R_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
allocate(R_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
allocate(R_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/read_mesh_databases.f90 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/read_mesh_databases.f90 2009-09-24 00:30:14 UTC (rev 15688)
@@ -52,6 +52,19 @@
read(27) rho_vp
read(27) rho_vs
read(27) iflag_attenuation_store
+
+ ! checks attenuation flags: see integers defined in constants.h
+ if( ATTENUATION ) then
+ if( minval(iflag_attenuation_store(:,:,:,:)) < 1 ) then
+ close(27)
+ call exit_MPI(myrank,'something is wrong with the mesh attenuation: flag entry is invalid')
+ endif
+ if( maxval(iflag_attenuation_store(:,:,:,:)) > NUM_REGIONS_ATTENUATION ) then
+ close(27)
+ call exit_MPI(myrank,'something is wrong with the mesh attenuation: flag entry exceeds number of defined attenuation flags in constants.h')
+ endif
+ endif
+
read(27) NSPEC2DMAX_XMIN_XMAX_ext
read(27) NSPEC2DMAX_YMIN_YMAX_ext
allocate(nimin(2,NSPEC2DMAX_YMIN_YMAX_ext),nimax(2,NSPEC2DMAX_YMIN_YMAX_ext),nkmin_eta(2,NSPEC2DMAX_YMIN_YMAX_ext))
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90 2009-09-24 00:30:14 UTC (rev 15688)
@@ -117,12 +117,16 @@
write(IOUT,*) '!'
! if attenuation is off, set dummy size of arrays to one
+! both parameters are obsolete for specfem3D
+! they are only used in ampuero_implicit_ABC_specfem3D.f90
+ write(IOUT,*) '! only needed for ampuero_implicit_ABC_specfem3D.f90 compilation: '
+ write(IOUT,*) '! (uncomment next line) '
if(ATTENUATION) then
-! write(IOUT,*) 'integer, parameter :: NSPEC_ATTENUATION = ', NSPEC_AB
- write(IOUT,*) 'logical, parameter :: ATTENUATION_VAL = .true.'
+ write(IOUT,*) '! integer, parameter :: NSPEC_ATTENUATION = ', NSPEC_AB
+! write(IOUT,*) '! logical, parameter :: ATTENUATION_VAL = .true.'
else
- write(IOUT,*) 'integer, parameter :: NSPEC_ATTENUATION = ', 1
- write(IOUT,*) 'logical, parameter :: ATTENUATION_VAL = .false.'
+ write(IOUT,*) '! integer, parameter :: NSPEC_ATTENUATION = ', 1
+! write(IOUT,*) '! logical, parameter :: ATTENUATION_VAL = .false.'
endif
write(IOUT,*)
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D_par.f90 2009-09-23 17:08:07 UTC (rev 15687)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D_par.f90 2009-09-24 00:30:14 UTC (rev 15688)
@@ -256,9 +256,12 @@
! timer MPI
double precision, external :: wtime
- integer ihours,iminutes,iseconds,int_tCPU
- double precision time_start,tCPU
+ integer :: ihours,iminutes,iseconds,int_tCPU, &
+ ihours_remain,iminutes_remain,iseconds_remain,int_t_remain, &
+ ihours_total,iminutes_total,iseconds_total,int_t_total
+ double precision :: time_start,tCPU,t_remain,t_total
+
! parameters read from parameter file
integer NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,SIMULATION_TYPE
integer NSOURCES
More information about the CIG-COMMITS
mailing list