[cig-commits] r17093 - in seismo/3D/SPECFEM3D/trunk: . EXAMPLES EXAMPLES/homogeneous_halfspace EXAMPLES/layered_halfspace EXAMPLES/tomographic_model EXAMPLES/waterlayered_halfspace decompose_mesh_SCOTCH
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Aug 16 13:37:29 PDT 2010
Author: danielpeter
Date: 2010-08-16 13:37:29 -0700 (Mon, 16 Aug 2010)
New Revision: 17093
Added:
seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/
seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/CMTSOLUTION
seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/Par_file
seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/README
seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/STATIONS
seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/boundary_definition.py
seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/create_tomography_model_file.sh
seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/cubit2specfem3d.py
seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/run_boundary_definition.py
seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/run_cubit2specfem3d.py
seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/tomoblock_mesh.py
Modified:
seismo/3D/SPECFEM3D/trunk/EXAMPLES/homogeneous_halfspace/cubit2specfem3d.py
seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py
seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/CMTSOLUTION
seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/Par_file
seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/README
seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/cubit2specfem3d.py
seismo/3D/SPECFEM3D/trunk/EXAMPLES/waterlayered_halfspace/cubit2specfem3d.py
seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90
seismo/3D/SPECFEM3D/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/decompose_mesh_SCOTCH/Makefile
seismo/3D/SPECFEM3D/trunk/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/generate_databases.f90
seismo/3D/SPECFEM3D/trunk/get_model.f90
seismo/3D/SPECFEM3D/trunk/initialize_simulation.f90
seismo/3D/SPECFEM3D/trunk/locate_source.f90
seismo/3D/SPECFEM3D/trunk/memory_eval.f90
seismo/3D/SPECFEM3D/trunk/model_tomography.f90
seismo/3D/SPECFEM3D/trunk/read_mesh_databases.f90
seismo/3D/SPECFEM3D/trunk/serial.f90
seismo/3D/SPECFEM3D/trunk/specfem3D_par.f90
seismo/3D/SPECFEM3D/trunk/write_VTK_data.f90
Log:
adds example for tomographic models; fixes bug in serial compilation (./configure -without-mpi)
Modified: seismo/3D/SPECFEM3D/trunk/EXAMPLES/homogeneous_halfspace/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/homogeneous_halfspace/cubit2specfem3d.py 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/homogeneous_halfspace/cubit2specfem3d.py 2010-08-16 20:37:29 UTC (rev 17093)
@@ -539,7 +539,7 @@
meshfile=open(mesh_name,'w')
print 'Writing '+mesh_name+'.....'
num_elems=cubit.get_hex_count()
- print num_elems
+ print ' number of elements:',str(num_elems)
meshfile.write(str(num_elems)+'\n')
num_write=0
for block,flag in zip(self.block_mat,self.block_flag):
@@ -554,7 +554,7 @@
txt=txt+('%10i %10i %10i %10i %10i %10i %10i %10i\n')% nodes[:]
meshfile.write(txt)
meshfile.close()
- print 'Ok',str(num_elems)
+ print 'Ok'
def material_write(self,mat_name):
mat=open(mat_name,'w')
print 'Writing '+mat_name+'.....'
@@ -569,6 +569,7 @@
print 'Writing '+nodecoord_name+'.....'
node_list=cubit.parse_cubit_list('node','all')
num_nodes=len(node_list)
+ print ' number of nodes:',str(num_nodes)
nodecoord.write('%10i\n' % num_nodes)
#
for node in node_list:
Modified: seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py 2010-08-16 20:37:29 UTC (rev 17093)
@@ -21,24 +21,32 @@
cubit.cmd('volume 1 move x 67000 y 67000 z -30000')
# create vertices for discontinuity
-cubit.cmd('split curve 9 distance 3000')
-cubit.cmd('split curve 10 distance 3000')
-cubit.cmd('split curve 11 distance 3000')
-cubit.cmd('split curve 12 distance 3000')
+distance = 3000
+cubit.cmd('split curve 9 distance '+str(distance))
+cubit.cmd('split curve 10 distance '+str(distance))
+cubit.cmd('split curve 11 distance '+str(distance))
+cubit.cmd('split curve 12 distance '+str(distance))
# create surface for interface
+# surface at 3 km depth
cubit.cmd('create surface vertex 9 10 12 11')
cubit.cmd('section volume 1 with surface 7 keep normal')
cubit.cmd('section volume 1 with surface 7 reverse')
# create vertices for auxiliary interface to allow for refinement
-cubit.cmd('split curve 29 distance 9000')
-cubit.cmd('split curve 31 distance 9000')
-cubit.cmd('split curve 32 distance 9000')
-cubit.cmd('split curve 36 distance 9000')
+#distance = 9000
+# to have a surface at 25050 m depth, such that point force source can be located exactly
+# on a GLL point at that depth as in Komatitsch et al. (1999)
+distance = 22050
+cubit.cmd('split curve 29 distance '+str(distance))
+cubit.cmd('split curve 31 distance '+str(distance))
+cubit.cmd('split curve 32 distance '+str(distance))
+cubit.cmd('split curve 36 distance '+str(distance))
+
# create surface for buffer interface to refine BELOW the discontinuity
+# surface at 3 km depth
cubit.cmd('create surface vertex 25 26 28 27')
cubit.cmd('section volume 3 with surface 19 keep normal')
@@ -50,15 +58,27 @@
cubit.cmd('imprint all')
# Meshing the volumes
-cubit.cmd('volume 3 size 3589.2')
+## middle volume
+#cubit.cmd('volume 3 size 3589.2')
+cubit.cmd('volume 3 size 4500.')
cubit.cmd('mesh volume 3')
+# refine boundary surface between top and middle volume
+# this will create a tripling layer from the bottom to the top in the middle volume
cubit.cmd('refine surface 8 numsplit 1 bias 1.0 depth 1')
-cubit.cmd('volume 1 size 1196.4')
+#cubit.cmd('pause')
+
+## top volume
+#cubit.cmd('volume 1 size 1196.4')
+cubit.cmd('volume 1 size 1000.0')
cubit.cmd('mesh volume 1')
-cubit.cmd('volume 5 size 4785.71')
+#cubit.cmd('pause')
+
+## bottom volume
+#cubit.cmd('volume 5 size 4785.71')
+cubit.cmd('volume 5 size 4500.')
cubit.cmd('mesh volume 5')
#### End of meshing
Modified: seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/CMTSOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/CMTSOLUTION 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/CMTSOLUTION 2010-08-16 20:37:29 UTC (rev 17093)
@@ -1,10 +1,10 @@
PDE 1999 01 01 00 00 00.00 67000 67000 -25000 4.2 4.2 FIG8
event name: FIG8_vertical
time shift: 0.0000
-half duration: 2.0
+half duration: 0.4
latitude: 67000.0
longitude: 67000.0
-depth: -25000.0
+depth: 25.05
Mrr: 1.000000e+23
Mtt: 0.000000
Mpp: 0.000000
Modified: seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/Par_file 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/Par_file 2010-08-16 20:37:29 UTC (rev 17093)
@@ -8,10 +8,10 @@
SUPPRESS_UTM_PROJECTION = .true.
# number of MPI processors
-NPROC = 4
+NPROC = 6
# time step parameters
-NSTEP = 4500
+NSTEP = 6000
DT = 0.0065d0
# parameters describing the model
@@ -22,7 +22,7 @@
ANISOTROPY = .false.
# absorbing boundary conditions for a regional simulation
-ABSORBING_CONDITIONS = .false.
+ABSORBING_CONDITIONS = .true.
# save AVS or OpenDX movies
MOVIE_SURFACE = .false.
Modified: seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/README
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/README 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/README 2010-08-16 20:37:29 UTC (rev 17093)
@@ -2,6 +2,8 @@
README
----------------------------------------------------------------------
+This example is used for validation with a layer-cake solution from Komatitsch et al. (1999),
+using a 2-layer model as shown in their Figure 8, left.
step-by-step tutorial:
@@ -70,13 +72,29 @@
4. run simulation:
+ - modify constants.h:
+ to use a vertical force source, with a Ricker wavelet source time function,
+ turn flag on for parameter USE_FORCE_POINT_SOURCE:
+ ...
+ logical, parameter :: USE_FORCE_POINT_SOURCE = .true.
+ ...
+
- compile specfem3D:
> make xspecfem3D
- submit job script:
> qsub go_solver_pbs.bash
+Optional:
+ Compare the solution traces:
+ 'OUTPUT_FILES/X31.BHZ.semd'
+ and
+ 'OUTPUT_FILES/X55.BHZ.semd'
+ with the solutions of Komatitsch et al. (1999, Figure 9, bottom and top) provided in:
+ 'EXAMPLES/VALIDATION_3D_SEM_SIMPLER_LAYER_SOURCE/REF_SEIS/Ux_file_Ycste_0031.txt'
+ and
+ 'EXAMPLES/VALIDATION_3D_SEM_SIMPLER_LAYER_SOURCE/REF_SEIS/Ux_file_Ycste_0055.txt'
\ No newline at end of file
Modified: seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/cubit2specfem3d.py 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/layered_halfspace/cubit2specfem3d.py 2010-08-16 20:37:29 UTC (rev 17093)
@@ -539,7 +539,7 @@
meshfile=open(mesh_name,'w')
print 'Writing '+mesh_name+'.....'
num_elems=cubit.get_hex_count()
- print num_elems
+ print ' number of elements:',str(num_elems)
meshfile.write(str(num_elems)+'\n')
num_write=0
for block,flag in zip(self.block_mat,self.block_flag):
@@ -554,7 +554,7 @@
txt=txt+('%10i %10i %10i %10i %10i %10i %10i %10i\n')% nodes[:]
meshfile.write(txt)
meshfile.close()
- print 'Ok',str(num_elems)
+ print 'Ok'
def material_write(self,mat_name):
mat=open(mat_name,'w')
print 'Writing '+mat_name+'.....'
@@ -569,6 +569,7 @@
print 'Writing '+nodecoord_name+'.....'
node_list=cubit.parse_cubit_list('node','all')
num_nodes=len(node_list)
+ print ' number of nodes:',str(num_nodes)
nodecoord.write('%10i\n' % num_nodes)
#
for node in node_list:
Added: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/CMTSOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/CMTSOLUTION (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/CMTSOLUTION 2010-08-16 20:37:29 UTC (rev 17093)
@@ -0,0 +1,13 @@
+PDE 1999 01 01 00 00 00.00 67000 67000 -25000 4.2 4.2 hom_explosion
+event name: hom_explosion
+time shift: 0.0000
+half duration: 5.0
+latitude: 67000.0
+longitude: 67000.0
+depth: 25.0
+Mrr: 1.000000e+23
+Mtt: 1.000000e+23
+Mpp: 1.000000e+23
+Mrt: 0.000000
+Mrp: 0.000000
+Mtp: 0.000000
Added: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/Par_file (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/Par_file 2010-08-16 20:37:29 UTC (rev 17093)
@@ -0,0 +1,51 @@
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint, 3 = both simultaneously
+SAVE_FORWARD = .false.
+
+# UTM projection parameters
+UTM_PROJECTION_ZONE = 11
+SUPPRESS_UTM_PROJECTION = .true.
+
+# number of MPI processors
+NPROC = 4
+
+# time step parameters
+NSTEP = 1000
+DT = 0.05d0
+
+# parameters describing the model
+OCEANS = .false.
+TOPOGRAPHY = .false.
+ATTENUATION = .false.
+USE_OLSEN_ATTENUATION = .false.
+ANISOTROPY = .false.
+
+# absorbing boundary conditions for a regional simulation
+ABSORBING_CONDITIONS = .false.
+
+# save AVS or OpenDX movies
+MOVIE_SURFACE = .false.
+MOVIE_VOLUME = .false.
+NTSTEP_BETWEEN_FRAMES = 200
+CREATE_SHAKEMAP = .false.
+SAVE_DISPLACEMENT = .false.
+USE_HIGHRES_FOR_MOVIES = .false.
+HDUR_MOVIE = 0.0
+
+# save AVS or OpenDX mesh files to check the mesh
+SAVE_MESH_FILES = .true.
+
+# path to store the local database file on each node
+LOCAL_PATH = DATABASES_MPI
+
+# interval at which we output time step info and max of norm of displacement
+NTSTEP_BETWEEN_OUTPUT_INFO = 500
+
+# interval in time steps for writing of seismograms
+NTSTEP_BETWEEN_OUTPUT_SEISMOS = 10000
+
+# print source time function
+PRINT_SOURCE_TIME_FUNCTION = .false.
+
+
Added: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/README
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/README (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/README 2010-08-16 20:37:29 UTC (rev 17093)
@@ -0,0 +1,119 @@
+----------------------------------------------------------------------
+README
+----------------------------------------------------------------------
+
+step-by-step tutorial:
+
+0. check that all software is available (or that modules are loaded):
+ intel/openmpi, cubit, scotch, python, gnuplot
+
+
+1. configure package:
+
+ - From the SPECFEM3D root directory SPECFEM3D/
+ configure the package, e.g. using intel's ifort compiler:
+ > cd SPECFEM3D
+ > ./configure F90=ifort
+
+ If successful, this will generate the files:
+ Makefile, constants.h, and precision.h, among others
+
+ - copy two run scripts from SPECFEM3D/UTILS/Cluster/
+ into SPECFEM3D/, e.g.,
+ pbs/go_generate_databases_pbs.bash
+ pbs/go_solver_pbs.bash
+
+
+2.a create mesh:
+
+ - change to the examples directory SPECFEM3D/EXAMPLES/tomographic_model:
+ > cd EXAMPLES/tomographic_model
+
+ - open the cubit GUI:
+ > claro (or cubit)
+
+ then run meshing script:
+ claro -> Menu "Tools" -> "Play Journal File" ... and select file: "tomoblock_mesh.py"
+
+ if everything goes fine, this creates all the mesh files in a subdirectory MESH/:
+ MESH/absorbing_surface_file_bottom
+ MESH/absorbing_surface_file_xmax
+ MESH/absorbing_surface_file_xmin
+ MESH/absorbing_surface_file_ymax
+ MESH/absorbing_surface_file_ymin
+ MESH/free_surface_file
+ MESH/materials_file
+ MESH/mesh_file
+ MESH/nodes_coords_file
+ MESH/nummaterial_velocity_file
+
+ the cubit graphics window should show a mesh similar to the file homogeneous.png
+
+2.b create tomography_model file:
+
+ in EXAMPLES/tomographic_model, run bash script:
+
+ > sh create_tomography_model_file.sh
+
+ which creates an example model: tomography_model.xyz
+ ( elastic model with a constant velocity gradient with depth)
+
+ copy this file into directory SPECFEM3D/DATA/
+
+
+3. decompose mesh files:
+
+ - compile decomposer in directory SPECFEM3D/decompose_mesh_SCOTCH/:
+ > make
+
+ NOTE 1: check that the two scotch libraries are properly specified in Makefile
+ NOTE 2: compile with the same compiler (ifort or gfortran) used
+ for the SCOTCH libraries
+
+ - run decomposer in directory
+ (example assumes 4 partitions with mesh files in ../EXAMPLES/homogeneous_halfspace/MESH/)
+ > ./xdecompose_mesh_SCOTCH 4 ../EXAMPLES/tomographic_model/MESH/ ../DATABASES_MPI/
+
+ this creates mesh partitions "proc000***_Database" in directory DATABASES_MPI/.
+ (you can then specify "DATABASES_MPI" in "Par_file" for your "LOCAL_PATH")
+
+
+4. generate databases:
+
+ - compile generate_databases from SPECFEM3D/ :
+ > cd SPECFEM3D
+ > make xgenerate_databases
+
+ - submit job script
+ > qsub go_generate_databases_pbs.bash
+
+ NOTE: this script will need to be tailored to your cluster, e.g.,
+ > bsub < go_generate_databases_lsf.bash
+
+ this will create binary mesh files, e.g. "proc000***_external_mesh.bin"
+ in directory DATABASES_MPI/.
+
+
+5. run simulation:
+
+ - copy three files -- Par_file CMTSOLUTION STATIONS -- from
+ SPECFEM3D/EXAMPLES/tomographic_model/ to SPECFEM3D/DATA/
+
+ - compile specfem3D:
+ > make xspecfem3D
+
+ - submit job script:
+ > qsub go_solver_pbs.bash
+
+ NOTE 1: this script will need to be tailored to your cluster, e.g.,
+ > bsub < go_solver_lsf.bash
+ NOTE 2: the simulation runs on 4 cores and should take about 15 minutes,
+ and you can track the progress with the timestamp files
+
+ - when the job is complete, you should have 3 sets (semd,semv,sema)
+ of 672 (ls -1 *semd | wc) in the directory OUTPUT_FILES,
+ as well as 3 timestamp****** files
+
+
+===========================================================
+
Added: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/STATIONS
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/STATIONS (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/STATIONS 2010-08-16 20:37:29 UTC (rev 17093)
@@ -0,0 +1,224 @@
+X1 DB 67000.00 0.000000 0.0 0.0
+X2 DB 67000.00 1196.429 0.0 0.0
+X3 DB 67000.00 2392.857 0.0 0.0
+X4 DB 67000.00 3589.286 0.0 0.0
+X5 DB 67000.00 4785.714 0.0 0.0
+X6 DB 67000.00 5982.143 0.0 0.0
+X7 DB 67000.00 7178.571 0.0 0.0
+X8 DB 67000.00 8375.000 0.0 0.0
+X9 DB 67000.00 9571.429 0.0 0.0
+X10 DB 67000.00 10767.86 0.0 0.0
+X11 DB 67000.00 11964.29 0.0 0.0
+X12 DB 67000.00 13160.71 0.0 0.0
+X13 DB 67000.00 14357.14 0.0 0.0
+X14 DB 67000.00 15553.57 0.0 0.0
+X15 DB 67000.00 16750.00 0.0 0.0
+X16 DB 67000.00 17946.43 0.0 0.0
+X17 DB 67000.00 19142.86 0.0 0.0
+X18 DB 67000.00 20339.29 0.0 0.0
+X19 DB 67000.00 21535.71 0.0 0.0
+X20 DB 67000.00 22732.14 0.0 0.0
+X21 DB 67000.00 23928.57 0.0 0.0
+X22 DB 67000.00 25125.00 0.0 0.0
+X23 DB 67000.00 26321.43 0.0 0.0
+X24 DB 67000.00 27517.86 0.0 0.0
+X25 DB 67000.00 28714.29 0.0 0.0
+X26 DB 67000.00 29910.71 0.0 0.0
+X27 DB 67000.00 31107.14 0.0 0.0
+X28 DB 67000.00 32303.57 0.0 0.0
+X29 DB 67000.00 33500.00 0.0 0.0
+X30 DB 67000.00 34696.43 0.0 0.0
+X31 DB 67000.00 35892.86 0.0 0.0
+X32 DB 67000.00 37089.29 0.0 0.0
+X33 DB 67000.00 38285.71 0.0 0.0
+X34 DB 67000.00 39482.14 0.0 0.0
+X35 DB 67000.00 40678.57 0.0 0.0
+X36 DB 67000.00 41875.00 0.0 0.0
+X37 DB 67000.00 43071.43 0.0 0.0
+X38 DB 67000.00 44267.86 0.0 0.0
+X39 DB 67000.00 45464.29 0.0 0.0
+X40 DB 67000.00 46660.71 0.0 0.0
+X41 DB 67000.00 47857.14 0.0 0.0
+X42 DB 67000.00 49053.57 0.0 0.0
+X43 DB 67000.00 50250.00 0.0 0.0
+X44 DB 67000.00 51446.43 0.0 0.0
+X45 DB 67000.00 52642.86 0.0 0.0
+X46 DB 67000.00 53839.29 0.0 0.0
+X47 DB 67000.00 55035.71 0.0 0.0
+X48 DB 67000.00 56232.14 0.0 0.0
+X49 DB 67000.00 57428.57 0.0 0.0
+X50 DB 67000.00 58625.00 0.0 0.0
+X51 DB 67000.00 59821.43 0.0 0.0
+X52 DB 67000.00 61017.86 0.0 0.0
+X53 DB 67000.00 62214.29 0.0 0.0
+X54 DB 67000.00 63410.71 0.0 0.0
+X55 DB 67000.00 64607.14 0.0 0.0
+X56 DB 67000.00 65803.57 0.0 0.0
+X57 DB 67000.00 67000.00 0.0 0.0
+X58 DB 67000.00 68196.43 0.0 0.0
+X59 DB 67000.00 69392.86 0.0 0.0
+X60 DB 67000.00 70589.29 0.0 0.0
+X61 DB 67000.00 71785.71 0.0 0.0
+X62 DB 67000.00 72982.14 0.0 0.0
+X63 DB 67000.00 74178.57 0.0 0.0
+X64 DB 67000.00 75375.00 0.0 0.0
+X65 DB 67000.00 76571.43 0.0 0.0
+X66 DB 67000.00 77767.86 0.0 0.0
+X67 DB 67000.00 78964.29 0.0 0.0
+X68 DB 67000.00 80160.71 0.0 0.0
+X69 DB 67000.00 81357.14 0.0 0.0
+X70 DB 67000.00 82553.57 0.0 0.0
+X71 DB 67000.00 83750.00 0.0 0.0
+X72 DB 67000.00 84946.43 0.0 0.0
+X73 DB 67000.00 86142.86 0.0 0.0
+X74 DB 67000.00 87339.29 0.0 0.0
+X75 DB 67000.00 88535.71 0.0 0.0
+X76 DB 67000.00 89732.14 0.0 0.0
+X77 DB 67000.00 90928.57 0.0 0.0
+X78 DB 67000.00 92125.00 0.0 0.0
+X79 DB 67000.00 93321.43 0.0 0.0
+X80 DB 67000.00 94517.86 0.0 0.0
+X81 DB 67000.00 95714.29 0.0 0.0
+X82 DB 67000.00 96910.71 0.0 0.0
+X83 DB 67000.00 98107.14 0.0 0.0
+X84 DB 67000.00 99303.57 0.0 0.0
+X85 DB 67000.00 100500.0 0.0 0.0
+X86 DB 67000.00 101696.4 0.0 0.0
+X87 DB 67000.00 102892.9 0.0 0.0
+X88 DB 67000.00 104089.3 0.0 0.0
+X89 DB 67000.00 105285.7 0.0 0.0
+X90 DB 67000.00 106482.1 0.0 0.0
+X91 DB 67000.00 107678.6 0.0 0.0
+X92 DB 67000.00 108875.0 0.0 0.0
+X93 DB 67000.00 110071.4 0.0 0.0
+X94 DB 67000.00 111267.9 0.0 0.0
+X95 DB 67000.00 112464.3 0.0 0.0
+X96 DB 67000.00 113660.7 0.0 0.0
+X97 DB 67000.00 114857.1 0.0 0.0
+X98 DB 67000.00 116053.6 0.0 0.0
+X99 DB 67000.00 117250.0 0.0 0.0
+X100 DB 67000.00 118446.4 0.0 0.0
+X101 DB 67000.00 119642.9 0.0 0.0
+X102 DB 67000.00 120839.3 0.0 0.0
+X103 DB 67000.00 122035.7 0.0 0.0
+X104 DB 67000.00 123232.1 0.0 0.0
+X105 DB 67000.00 124428.6 0.0 0.0
+X106 DB 67000.00 125625.0 0.0 0.0
+X107 DB 67000.00 126821.4 0.0 0.0
+X108 DB 67000.00 128017.9 0.0 0.0
+X109 DB 67000.00 129214.3 0.0 0.0
+X110 DB 67000.00 130410.7 0.0 0.0
+X111 DB 67000.00 131607.1 0.0 0.0
+X112 DB 67000.00 132803.6 0.0 0.0
+Y1 DB 0.000000 67000.00 0.0 0.0
+Y2 DB 1196.429 67000.00 0.0 0.0
+Y3 DB 2392.857 67000.00 0.0 0.0
+Y4 DB 3589.286 67000.00 0.0 0.0
+Y5 DB 4785.714 67000.00 0.0 0.0
+Y6 DB 5982.143 67000.00 0.0 0.0
+Y7 DB 7178.571 67000.00 0.0 0.0
+Y8 DB 8375.000 67000.00 0.0 0.0
+Y9 DB 9571.429 67000.00 0.0 0.0
+Y10 DB 10767.86 67000.00 0.0 0.0
+Y11 DB 11964.29 67000.00 0.0 0.0
+Y12 DB 13160.71 67000.00 0.0 0.0
+Y13 DB 14357.14 67000.00 0.0 0.0
+Y14 DB 15553.57 67000.00 0.0 0.0
+Y15 DB 16750.00 67000.00 0.0 0.0
+Y16 DB 17946.43 67000.00 0.0 0.0
+Y17 DB 19142.86 67000.00 0.0 0.0
+Y18 DB 20339.29 67000.00 0.0 0.0
+Y19 DB 21535.71 67000.00 0.0 0.0
+Y20 DB 22732.14 67000.00 0.0 0.0
+Y21 DB 23928.57 67000.00 0.0 0.0
+Y22 DB 25125.00 67000.00 0.0 0.0
+Y23 DB 26321.43 67000.00 0.0 0.0
+Y24 DB 27517.86 67000.00 0.0 0.0
+Y25 DB 28714.29 67000.00 0.0 0.0
+Y26 DB 29910.71 67000.00 0.0 0.0
+Y27 DB 31107.14 67000.00 0.0 0.0
+Y28 DB 32303.57 67000.00 0.0 0.0
+Y29 DB 33500.00 67000.00 0.0 0.0
+Y30 DB 34696.43 67000.00 0.0 0.0
+Y31 DB 35892.86 67000.00 0.0 0.0
+Y32 DB 37089.29 67000.00 0.0 0.0
+Y33 DB 38285.71 67000.00 0.0 0.0
+Y34 DB 39482.14 67000.00 0.0 0.0
+Y35 DB 40678.57 67000.00 0.0 0.0
+Y36 DB 41875.00 67000.00 0.0 0.0
+Y37 DB 43071.43 67000.00 0.0 0.0
+Y38 DB 44267.86 67000.00 0.0 0.0
+Y39 DB 45464.29 67000.00 0.0 0.0
+Y40 DB 46660.71 67000.00 0.0 0.0
+Y41 DB 47857.14 67000.00 0.0 0.0
+Y42 DB 49053.57 67000.00 0.0 0.0
+Y43 DB 50250.00 67000.00 0.0 0.0
+Y44 DB 51446.43 67000.00 0.0 0.0
+Y45 DB 52642.86 67000.00 0.0 0.0
+Y46 DB 53839.29 67000.00 0.0 0.0
+Y47 DB 55035.71 67000.00 0.0 0.0
+Y48 DB 56232.14 67000.00 0.0 0.0
+Y49 DB 57428.57 67000.00 0.0 0.0
+Y50 DB 58625.00 67000.00 0.0 0.0
+Y51 DB 59821.43 67000.00 0.0 0.0
+Y52 DB 61017.86 67000.00 0.0 0.0
+Y53 DB 62214.29 67000.00 0.0 0.0
+Y54 DB 63410.71 67000.00 0.0 0.0
+Y55 DB 64607.14 67000.00 0.0 0.0
+Y56 DB 65803.57 67000.00 0.0 0.0
+Y57 DB 67000.00 67000.00 0.0 0.0
+Y58 DB 68196.43 67000.00 0.0 0.0
+Y59 DB 69392.86 67000.00 0.0 0.0
+Y60 DB 70589.29 67000.00 0.0 0.0
+Y61 DB 71785.71 67000.00 0.0 0.0
+Y62 DB 72982.14 67000.00 0.0 0.0
+Y63 DB 74178.57 67000.00 0.0 0.0
+Y64 DB 75375.00 67000.00 0.0 0.0
+Y65 DB 76571.43 67000.00 0.0 0.0
+Y66 DB 77767.86 67000.00 0.0 0.0
+Y67 DB 78964.29 67000.00 0.0 0.0
+Y68 DB 80160.71 67000.00 0.0 0.0
+Y69 DB 81357.14 67000.00 0.0 0.0
+Y70 DB 82553.57 67000.00 0.0 0.0
+Y71 DB 83750.00 67000.00 0.0 0.0
+Y72 DB 84946.43 67000.00 0.0 0.0
+Y73 DB 86142.86 67000.00 0.0 0.0
+Y74 DB 87339.29 67000.00 0.0 0.0
+Y75 DB 88535.71 67000.00 0.0 0.0
+Y76 DB 89732.14 67000.00 0.0 0.0
+Y77 DB 90928.57 67000.00 0.0 0.0
+Y78 DB 92125.00 67000.00 0.0 0.0
+Y79 DB 93321.43 67000.00 0.0 0.0
+Y80 DB 94517.86 67000.00 0.0 0.0
+Y81 DB 95714.29 67000.00 0.0 0.0
+Y82 DB 96910.71 67000.00 0.0 0.0
+Y83 DB 98107.14 67000.00 0.0 0.0
+Y84 DB 99303.57 67000.00 0.0 0.0
+Y85 DB 100500.0 67000.00 0.0 0.0
+Y86 DB 101696.4 67000.00 0.0 0.0
+Y87 DB 102892.9 67000.00 0.0 0.0
+Y88 DB 104089.3 67000.00 0.0 0.0
+Y89 DB 105285.7 67000.00 0.0 0.0
+Y90 DB 106482.1 67000.00 0.0 0.0
+Y91 DB 107678.6 67000.00 0.0 0.0
+Y92 DB 108875.0 67000.00 0.0 0.0
+Y93 DB 110071.4 67000.00 0.0 0.0
+Y94 DB 111267.9 67000.00 0.0 0.0
+Y95 DB 112464.3 67000.00 0.0 0.0
+Y96 DB 113660.7 67000.00 0.0 0.0
+Y97 DB 114857.1 67000.00 0.0 0.0
+Y98 DB 116053.6 67000.00 0.0 0.0
+Y99 DB 117250.0 67000.00 0.0 0.0
+Y100 DB 118446.4 67000.00 0.0 0.0
+Y101 DB 119642.9 67000.00 0.0 0.0
+Y102 DB 120839.3 67000.00 0.0 0.0
+Y103 DB 122035.7 67000.00 0.0 0.0
+Y104 DB 123232.1 67000.00 0.0 0.0
+Y105 DB 124428.6 67000.00 0.0 0.0
+Y106 DB 125625.0 67000.00 0.0 0.0
+Y107 DB 126821.4 67000.00 0.0 0.0
+Y108 DB 128017.9 67000.00 0.0 0.0
+Y109 DB 129214.3 67000.00 0.0 0.0
+Y110 DB 130410.7 67000.00 0.0 0.0
+Y111 DB 131607.1 67000.00 0.0 0.0
+Y112 DB 132803.6 67000.00 0.0 0.0
Added: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/boundary_definition.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/boundary_definition.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/boundary_definition.py 2010-08-16 20:37:29 UTC (rev 17093)
@@ -0,0 +1,378 @@
+#############################################################################
+# boundary_definition.py #
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+
+def define_absorbing_surf():
+ """
+ define the absorbing surfaces for a layered topological box where boundary are surfaces parallel to the axis.
+ it returns absorbing_surf,absorbing_surf_xmin,absorbing_surf_xmax,absorbing_surf_ymin,absorbing_surf_ymax,absorbing_surf_bottom,topo_surf
+ where
+ absorbing_surf is the list of all the absorbing boundary surf
+ absorbing_surf_xmin is the list of the absorbing boundary surfaces that correnspond to x=xmin
+ ...
+ absorbing_surf_bottom is the list of the absorbing boundary surfaces that correspond to z=zmin
+ """
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ absorbing_surf=[]
+ absorbing_surf_xmin=[]
+ absorbing_surf_xmax=[]
+ absorbing_surf_ymin=[]
+ 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]
+ zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6] #it is the z_min of the box ... box= xmin,xmax,d,ymin,ymax,d,zmin...
+ xmin_box=cubit.get_total_bounding_box("volume",list_vol)[0]
+ xmax_box=cubit.get_total_bounding_box("volume",list_vol)[1]
+ 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)/x_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_xmin.append(k)
+ absorbing_surf.append(k)
+ 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)/y_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_ymin.append(k)
+ absorbing_surf.append(k)
+ 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)/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)/z_len)
+ normal=cubit.get_surface_normal(k)
+ zn=normal[2]
+ dn=abs(zn-1)
+ 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():
+ """
+ define the absorbing surfaces for a layered topological box where boundary surfaces are not parallel to the axis.
+ it returns absorbing_surf,topo_surf
+ where
+ absorbing_surf is the list of all the absorbing boundary surf
+ """
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ from sets import Set
+ def product(*args, **kwds):
+ # product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
+ # product(range(2), repeat=3) --> 000 001 010 011 100 101 110 111
+ pools = map(tuple, args) * kwds.get('repeat', 1)
+ result = [[]]
+ for pool in pools:
+ result = [x+[y] for x in result for y in pool]
+ return result
+ absorbing_surf=[]
+ absorbing_surf_xmin=[]
+ absorbing_surf_xmax=[]
+ absorbing_surf_ymin=[]
+ 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]
+ zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6] #it is the z_min of the box ... box= xmin,xmax,d,ymin,ymax,d,zmin...
+ xmin_box=cubit.get_total_bounding_box("volume",list_vol)[0]
+ xmax_box=cubit.get_total_bounding_box("volume",list_vol)[1]
+ 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")
+ lv=[]
+ for k in list_surf:
+ sbox=cubit.get_bounding_box('surface',k)
+ dzmax=abs((sbox[7] - zmax_box)/zmax_box)
+ dzmin=abs((sbox[6] - zmin_box)/zmin_box)
+ normal=cubit.get_surface_normal(k)
+ zn=normal[2]
+ if dzmax <= 0.001 and zn > 0.7:
+ top_surf.append(k)
+ list_vertex=cubit.get_relatives('surface',k,'vertex')
+ for v in list_vertex:
+ valence=cubit.get_valence(v)
+ if valence <= 4: #valence 3 is a corner, 4 is a vertex between 2 volumes, > 4 is a vertex not in the boundaries
+ lv.append(v)
+ elif dzmin <= 0.001 and zn < -0.7:
+ absorbing_surf.append(k)
+ lp=[]
+ combs=product(lv,lv)
+ for comb in combs:
+ v1=comb[0]
+ v2=comb[1]
+ c=Set(cubit.get_relatives("vertex",v1,"curve")) & Set(cubit.get_relatives("vertex",v2,"curve"))
+ if len(c) == 1:
+ p=cubit.get_center_point("curve",list(c)[0])
+ lp.append(p)
+ for k in list_surf:
+ center_point = cubit.get_center_point("surface", k)
+ for p in lp:
+ if abs((center_point[0] - p[0])/p[0]) <= 0.005 and abs((center_point[1] - p[1])/p[1]) <= 0.005:
+ absorbing_surf.append(k)
+ break
+ return absorbing_surf,top_surf
+
+def define_absorbing_surf_sphere():
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ surf=[]
+ list_surf=cubit.parse_cubit_list("surface","all")
+ for s in list_surf:
+ v=cubit.get_relatives('surface',s,'volume')
+ if len(v) == 1:
+ surf.append(s)
+ return surf
+
+def define_block():
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ list_vol=cubit.parse_cubit_list("volume","all")
+ init_n_vol=len(list_vol)
+ list_name=map(lambda x: 'vol'+x,map(str,list_vol))
+ return list_vol,list_name
+
+def build_block(vol_list,name):
+ from sets import Set
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ block_list=cubit.get_block_id_list()
+ if len(block_list) > 0:
+ id_block=max(block_list)
+ else:
+ id_block=0
+ for v,n in zip(vol_list,name):
+ id_block+=1
+ v_other=Set(vol_list)-Set([v])
+ #command= 'block '+str(id_block)+' hex in node in vol '+str(v)+' except hex in vol '+str(list(v_other))
+ command= 'block '+str(id_block)+' hex in vol '+str(v)
+ command = command.replace("["," ").replace("]"," ")
+ cubit.cmd(command)
+ command = "block "+str(id_block)+" name '"+n+"'"
+ cubit.cmd(command)
+
+def build_block_side(surf_list,name,obj='surface'):
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ 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))
+ txt2="block "+str(id_block)+" name '"+name+"'"
+ txt1=txt1.replace("["," ").replace("]"," ")
+ elif obj == 'node':
+ txt=obj+' in surface'
+ txt1= 'nodeset '+str(id_nodeset)+ ' '+ txt +' '+str(list(surf_list))
+ txt1 = txt1.replace("["," ").replace("]"," ")
+ txt2 = "nodeset "+str(id_nodeset)+" name '"+name+"'"
+ elif obj == 'face' or obj == 'edge':
+ txt=obj+' in surface'
+ txt1= 'block '+str(id_block)+ ' '+ txt +' '+str(list(surf_list))
+ txt1 = txt1.replace("["," ").replace("]"," ")
+ txt2 = "block "+str(id_block)+" name '"+name+"'"
+ else:
+ txt1=''
+ # 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)
+
+def define_bc(*args,**keys):
+ 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:
+ surf,topo=define_absorbing_surf_nopar()
+ v_list,name_list=define_block()
+ build_block(v_list,name_list)
+ 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:
+ build_block_side(xmin,entity+'_abs_xmin',obj=entity)
+ build_block_side(xmax,entity+'_abs_xmax',obj=entity)
+ build_block_side(ymin,entity+'_abs_ymin',obj=entity)
+ 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)
+ entities=args[0]
+ for entity in entities:
+ build_block_side(surf,entity+'_closedvol',obj=entity)
+
+
+
+#entities=['surface','face','edge','node','hex']
+#define_bc(entities,parallel=True)
+#define_bc(entities,parallel=False)
+#define_bc(entities,parallel=False,closed=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'
+#block 2 attribute count 3
+#block 3 attribute count 5
+#block 1 attribute index 2 1500
+#block 1 attribute count 2
+#block 3 attribute index 1 2
+#block 3 attribute index 2 5800
+#block 3 attribute index 3 3900
+#block 3 attribute index 4 1500
+#block 3 attribute index 5 3.5
+#block 1 name 'top'
+#block 3 name 'bottom'
+#block 2 attribute index 1 -1
+#block 2 attribute index 2 2
+#block 2 attribute count 4
+#block 2 attribute index 2 1
+#block 2 attribute index 3 2
+
+
+
+
Added: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/create_tomography_model_file.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/create_tomography_model_file.sh (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/create_tomography_model_file.sh 2010-08-16 20:37:29 UTC (rev 17093)
@@ -0,0 +1,55 @@
+#!/bin/sh
+
+# creates a simple, formatted tomography model with a constant velocity gradient
+# for a block model with dimensions 134000 x 134000 x 60000
+
+# origin points
+ORIG_X=0.
+ORIG_Y=0.
+ORIG_Z=0.
+
+# end points
+END_X=134000.
+END_Y=134000.
+END_Z=-60000. # depth in negative z-direction
+
+# spacing of given tomography points
+SPACING_X=2000.
+SPACING_Y=2000.
+SPACING_Z=-2000.
+
+# number of cell increments
+NX=68
+NY=68
+NZ=31
+
+# min/max values
+VP_MIN=2500.
+VP_MAX=8500.
+VS_MIN=1500.
+VS_MAX=7500.
+RHO_MIN=1500.
+RHO_MAX=1500.
+
+
+# header info
+echo "creating header info..."
+
+echo "$ORIG_X $ORIG_Y $ORIG_Z $END_X $END_Y $END_Z " > tmp.xyz
+echo "$SPACING_X $SPACING_Y $SPACING_Z " >> tmp.xyz
+echo "$NX $NY $NZ " >> tmp.xyz
+echo "$VP_MIN $VP_MAX $VS_MIN $VS_MAX $RHO_MIN $RHO_MAX" >> tmp.xyz
+
+# velocity gradient
+GRADIENT=0.1
+
+# adds point location and velocity model values
+echo "adding model values..."
+
+# format: lists first all x, then y, then z
+echo "1" | awk '{ for(k=0;k<NZ;k++){ for(j=0;j<NY;j++){for(i=0;i<NX;i++){ x=i*SPACING_X;y=j*SPACING_Y;z=k*SPACING_Z;vp=VP_MIN+GRADIENT*(-z);vs=VS_MIN + GRADIENT*(-z); rho=RHO_MIN;print x,y,z,vp,vs,rho }}} }' \
+ NX=$NX NY=$NY NZ=$NZ SPACING_X=$SPACING_X SPACING_Y=$SPACING_Y SPACING_Z=$SPACING_Z VP_MIN=$VP_MIN VS_MIN=$VS_MIN RHO_MIN=$RHO_MIN GRADIENT=$GRADIENT >> tmp.xyz
+
+mv tmp.xyz tomography_model.xyz
+
+echo "created file: tomography_model.xyz"
Property changes on: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/create_tomography_model_file.sh
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/cubit2specfem3d.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/cubit2specfem3d.py 2010-08-16 20:37:29 UTC (rev 17093)
@@ -0,0 +1,776 @@
+#!python
+#############################################################################
+# cubit2specfem3d.py #
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+#
+#for a complete definition of the format of the mesh in SPECFEM3D_SESAME check the manual (http:/):
+#
+#USAGE
+#
+#############################################################################
+#PREREQUISITE
+#The mesh must be prepared
+# automatically using the module boundary_definition (see boundary_definition.py for more information)
+#or
+# manually following the convention:
+# - each material should have a block defined by:
+# material domain_flag (acoustic/elastic/poroelastic)name,flag of the material (integer),p velocity
+# (or the full description: name, flag, vp, vs, rho, Q ... if not present these last 3 parameters will be
+# interpolated by module mat_parameter)
+# - each mesh should have the block definition for the face on the free_surface (topography),
+# the name of this block must be 'face_topo' or you can change the default name in mesh.topo defined in profile.
+# - each mesh should have the block definition for the faces on the absorbing boundaries,
+# one block for each surface with x=Xmin,x=Xmax,y=Ymin,y=Ymax and z=bottom. The names of
+# the blocks should contain the strings "xmin,xmax,ymin,ymax,bottom"
+#
+#############################################################################
+#RUN
+#In a python script or in the cubit python tab call:
+#
+# export2SESAME(path_exporting_mesh_SPECFEM3D_SESAME)
+#
+#the modele create a python class for the mesh: ex. profile=mesh()
+#and it export the files of the mesh needed by the partitioner of SESAME
+#
+#############################################################################
+#OUTPUT
+#The default output are 11 ASCII files:
+#__________________________________________________________________________________________
+#mesh_name='mesh_file' -> the file that contains the connectity of the all mesh
+# format:
+# number of elements
+# id_elements id_node1 id_node2 id_node3 id_node4 id_node5 id_node6 id_node7 id_node8
+# .....
+#
+#__________________________________________________________________________________________
+##nodecoord_name='nodes_coords_file' -> the file that contains the coordinates of the nodes of the all mesh
+# format:
+# number of nodes
+# id_node x_coordinate y_coordinate z_coordinate
+# .....
+#
+#__________________________________________________________________________________________
+##material_name='materials_file' -> the file that contains the material flag of the elements
+# format:
+# id_element flag
+# .....
+#
+#__________________________________________________________________________________________
+##nummaterial_name='nummaterial_velocity_file' -> table of the material properties
+# format:
+# flag rho vp vs 0 0 #full definition of the properties, flag > 0
+# .....
+# flag 'tomography' file_name #for interpolation with tomography
+# .....
+# flag 'interface' file_name flag_for_the_gll_below_the_interface
+# flag_for_the_gll_above_the_interface #for interpolation with interface
+#__________________________________________________________________________________________
+##absname='absorbing_surface_file' -> this file contains all the face in all the absorbing boundaries
+##absname_local='absorbing_surface_file'+'_xmin' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmin
+##absname_local='absorbing_surface_file'+'_xmax' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmax
+##absname_local='absorbing_surface_file'+'_ymin' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymin
+##absname_local='absorbing_surface_file'+'_ymax' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymax
+##absname_local='absorbing_surface_file'+'_bottom' -> this file contains all the face in the
+# absorbing boundary defined by z=bottom
+# format:
+# number of faces
+# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
+# ....
+#
+#__________________________________________________________________________________________
+##freename='free_surface_file' -> file with the hex on the free surface (usually the topography)
+# format:
+# number of faces
+# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
+#
+#__________________________________________________________________________________________
+# it is possible save only one (or more) file singularly: for example if you want only the nodecoord_file
+# call the module mesh.nodescoord_write(full path name)
+#
+#############################################################################
+
+import cubit
+
+class mtools(object):
+ """docstring for ciao"""
+ def __init__(self,frequency,list_surf,list_vp):
+ super(mtools, self).__init__()
+ self.frequency = frequency
+ self.list_surf = list_surf
+ self.list_vp = list_vp
+ self.ngll=5
+ self.percent_gll=0.172
+ self.point_wavelength=5
+ def __repr__(self):
+ txt='Meshing for frequency up to '+str(self.frequency)+'Hz\n'
+ for surf,vp in zip(self.list_surf,self.list_vp):
+ txt=txt+'surface '+str(surf)+', vp ='+str(vp)+' -> size '+str(self.freq2meshsize(vp)[0])\
+ +' -> dt '+str(self.freq2meshsize(vp)[0])+'\n'
+ return txt
+ def freq2meshsize(self,vp):
+ velocity=vp*.5
+ self.size=(1/2.5)*velocity/self.frequency*(self.ngll-1)/self.point_wavelength
+ self.dt=.4*self.size/vp*self.percent_gll
+ return self.size,self.dt
+ def mesh_it(self):
+ for surf,vp in zip(self.list_surf,self.list_vp):
+ command = "surface "+str(surf)+" size "+str(self.freq2meshsize(vp)[0])
+ cubit.cmd(command)
+ command = "surface "+str(surf)+ 'scheme pave'
+ cubit.cmd(command)
+ command = "mesh surf "+str(surf)
+ cubit.cmd(command)
+
+class block_tools:
+ def __int__(self):
+ pass
+ def create_blocks(self,mesh_entity,list_entity=None,):
+ if mesh_entity =='surface':
+ txt=' face in surface '
+ elif mesh_entity == 'curve':
+ txt=' edge in curve '
+ elif mesh_entity == 'group':
+ txt=' face in group '
+ if list_entity:
+ if not isinstance(list_entity,list):
+ list_entity=[list_entity]
+ for entity in list_entity:
+ iblock=cubit.get_next_block_id()
+ command = "block "+str(iblock)+ txt +str(entity)
+ cubit.cmd(command)
+ def material_file(self,filename):
+ matfile=open(filename,'w')
+ material=[]
+ for record in matfile:
+ mat_name,vp_str=record.split()
+ vp=float(vp_str)
+ material.append([mat_name,vp])
+ self.material=dict(material)
+ def assign_block_material(self,id_block,mat_name,vp=None):
+ try:
+ material=self.material
+ except:
+ material=None
+ cubit.cmd('block '+str(id_block)+' attribute count 2')
+ cubit.cmd('block '+str(id_block)+' attribute index 1 '+str(id_block))
+ if material:
+ if material.has_key(mat_name):
+ cubit.cmd('block '+str(id_block)+' attribute index 2 '+str(material[mat_name]))
+ print 'block '+str(id_block)+' - material '+mat_name+' - vp '+str(material[mat_name])+' from database'
+ elif vp:
+ cubit.cmd('block '+str(id_block)+' attribute index 2 '+str(vp))
+ print 'block '+str(id_block)+' - material '+mat_name+' - vp '+str(vp)
+ else:
+ print 'assignment impossible: check if '+mat_name+' is in the database or specify vp'
+
+class mesh_tools(block_tools):
+ """Tools for the mesh
+ #########
+ dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
+ #########
+ length=edge_length(edge)
+ return the length of a edge
+ #########
+ edge_min,length=edge_min_length(surface)
+ given the cubit id of a surface, it return the edge with minimun length
+ #########
+ """
+ def __int__(self):
+ pass
+ def seismic_resolution(self,edges,velocity,bins_d=None,bins_u=None,sidelist=None):
+ """
+ dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
+ """
+ ratiostore=1e10
+ dtstore=1e10
+ edgedtstore=-1
+ edgeratiostore=-1
+ for edge in edges:
+ d=self.edge_length(edge)
+ ratio=(1/2.5)*velocity/d*(self.ngll-1)/self.point_wavelength
+ dt=.4*d/velocity*self.percent_gll
+ if dt<dtstore:
+ dtstore=dt
+ edgedtstore=edge
+ if ratio < ratiostore:
+ ratiostore=ratio
+ edgeratiostore=edge
+ try:
+ for bin_d,bin_u,side in zip(bins_d,bins_u,sidelist):
+ if ratio >= bin_d and ratio < bin_u:
+ command = "sideset "+str(side)+" edge "+str(edge)
+ cubit.cmd(command)
+ #print command
+ break
+ except:
+ pass
+ return dtstore,edgedtstore,ratiostore,edgeratiostore
+ def edge_length(self,edge):
+ """
+ length=edge_length(edge)
+ return the length of a edge
+ """
+ from math import sqrt
+ nodes=cubit.get_connectivity('Edge',edge)
+ x0,y0,z0=cubit.get_nodal_coordinates(nodes[0])
+ x1,y1,z1=cubit.get_nodal_coordinates(nodes[1])
+ d=sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
+ return d
+ def edge_min_length(self,surface):
+ """
+ edge_min,length=edge_min_length(surface)
+ given the cubit id of a surface, it return the edge with minimun length
+ """
+ from math import sqrt
+ self.dmin=99999
+ edge_store=0
+ command = "group 'list_edge' add edge in surf "+str(surface)
+ command = command.replace("["," ").replace("]"," ")
+ #print command
+ cubit.cmd(command)
+ group=cubit.get_id_from_name("list_edge")
+ edges=cubit.get_group_edges(group)
+ command = "delete group "+ str(group)
+ cubit.cmd(command)
+ for edge in edges:
+ d=self.edge_length(edge)
+ if d<dmin:
+ self.dmin=d
+ edge_store=edge
+ self.edgemin=edge_store
+ return self.edgemin,self.dmin
+ def normal_check(self,nodes,normal):
+ tres=.2
+ p0=cubit.get_nodal_coordinates(nodes[0])
+ p1=cubit.get_nodal_coordinates(nodes[1])
+ p2=cubit.get_nodal_coordinates(nodes[2])
+ a=[p1[0]-p0[0],p1[1]-p0[1],p1[2]-p0[2]]
+ b=[p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]]
+ axb=[a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]]
+ dot=0.0
+ for i in (0,1,2):
+ dot=dot+axb[i]*normal[i]
+ if dot > 0:
+ return nodes
+ elif dot < 0:
+ return nodes[0],nodes[3],nodes[2],nodes[1]
+ else:
+ print 'error: surface normal, dot=0', axb,normal,dot,p0,p1,p2
+ def mesh_analysis(self,frequency):
+ from sets import Set
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ bins_d=[0.0001]+range(0,int(frequency)+1)+[1000]
+ bins_u=bins_d[1:]
+ dt=[]
+ ed_dt=[]
+ r=[]
+ ed_r=[]
+ nstart=cubit.get_next_sideset_id()
+ command = "del sideset all"
+ cubit.cmd(command)
+ for bin_d,bin_u in zip(bins_d,bins_u):
+ nsideset=cubit.get_next_sideset_id()
+ command='create sideset '+str(nsideset)
+ cubit.cmd(command)
+ command = "sideset "+str(nsideset)+ " name "+ "'ratio-["+str(bin_d)+"_"+str(bin_u)+"['"
+ cubit.cmd(command)
+ nend=cubit.get_next_sideset_id()
+ sidelist=range(nstart,nend)
+ for block in self.block_mat:
+ name=cubit.get_exodus_entity_name('block',block)
+ velocity=self.material[name][1]
+ if velocity > 0:
+ faces=cubit.get_block_faces(block)
+ edges=[]
+ for face in faces:
+ es=cubit.get_sub_elements("face", face, 1)
+ edges=edges+list(es)
+ edges=Set(edges)
+ dtstore,edgedtstore,ratiostore,edgeratiostore=self.seismic_resolution(edges,\
+ velocity,bins_d,bins_u,sidelist)
+ dt.append(dtstore)
+ ed_dt.append(edgedtstore)
+ r.append(ratiostore)
+ ed_r.append(edgeratiostore)
+ self.ddt=zip(ed_dt,dt)
+ self.dr=zip(ed_r,r)
+ def sorter(x, y):
+ return cmp(x[1],y[1])
+ self.ddt.sort(sorter)
+ self.dr.sort(sorter)
+ print self.ddt,self.dr
+ print 'Deltat minimum => edge:'+str(self.ddt[0][0])+' dt: '+str(self.ddt[0][1])
+ print 'minimum frequency resolved => edge:'+str(self.dr[0][0])+' frequency: '+str(self.dr[0][1])
+ return self.ddt[0],self.dr[0]
+
+class mesh(object,mesh_tools):
+ def __init__(self):
+ super(mesh, self).__init__()
+ self.mesh_name='mesh_file'
+ self.nodecoord_name='nodes_coords_file'
+ self.material_name='materials_file'
+ self.nummaterial_name='nummaterial_velocity_file'
+ self.absname='absorbing_surface_file'
+ self.freename='free_surface_file'
+ self.recname='STATIONS'
+ self.face='QUAD4'
+ self.face2='SHELL4'
+ self.hex='HEX8'
+ self.edge='BAR2'
+ self.topo='face_topo'
+ self.rec='receivers'
+ self.ngll=5
+ self.percent_gll=0.172
+ self.point_wavelength=5
+ self.block_definition()
+ cubit.cmd('compress')
+ def __repr__(self):
+ pass
+ def block_definition(self):
+ block_flag=[]
+ block_mat=[]
+ block_bc=[]
+ block_bc_flag=[]
+ material={}
+ bc={}
+ blocks=cubit.get_block_id_list()
+ for block in blocks:
+ name=cubit.get_exodus_entity_name('block',block)
+ type=cubit.get_block_element_type(block)
+ print block,name,blocks,type,self.hex,self.face
+ # block has hexahedral elements (HEX8)
+ if type == self.hex:
+ flag=None
+ vel=None
+ vs=None
+ rho=None
+ q=0
+ ani=0
+ # material domain id
+ if name.find("acoustic") >= 0 :
+ imaterial = 1
+ elif name.find("elastic") >= 0 :
+ imaterial = 2
+ elif name.find("poroelastic") >= 0 :
+ imaterial = 3
+ else :
+ imaterial = 0
+ print "block: ",name
+ print " could not find appropriate material for this block..."
+ print ""
+ break
+
+ nattrib=cubit.get_block_attribute_count(block)
+ if nattrib != 0:
+ # material flag:
+ # positive => material properties,
+ # negative => interface/tomography domain
+ flag=int(cubit.get_block_attribute_value(block,0))
+ if flag > 0 and nattrib >= 2:
+ # vp
+ vel=cubit.get_block_attribute_value(block,1)
+ if nattrib >= 3:
+ # vs
+ vs=cubit.get_block_attribute_value(block,2)
+ if nattrib >= 4:
+ #density
+ rho=cubit.get_block_attribute_value(block,3)
+ if nattrib >= 5:
+ #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
+ if nattrib == 6:
+ #anisotropy_flag
+ ani=cubit.get_block_attribute_value(block,5)
+ elif flag < 0:
+ # velocity model
+ vel=name
+ attrib=cubit.get_block_attribute_value(block,1)
+ if attrib == 1:
+ kind='interface'
+ flag_down=cubit.get_block_attribute_value(block,2)
+ flag_up=cubit.get_block_attribute_value(block,3)
+ elif attrib == 2:
+ kind='tomography'
+ else:
+ flag=block
+ vel,vs,rho,q,ani=(name,0,0,0,0)
+ block_flag.append(int(flag))
+ block_mat.append(block)
+ if flag > 0:
+ par=tuple([imaterial,flag,vel,vs,rho,q,ani])
+ elif flag < 0:
+ if kind=='interface':
+ par=tuple([imaterial,flag,kind,name,flag_down,flag_up])
+ elif kind=='tomography':
+ par=tuple([imaterial,flag,kind,name])
+ elif flag==0:
+ par=tuple([imaterial,flag,name])
+ material[block]=par
+ elif (type == self.face) or (type == self.face2) :
+ # block has surface elements (QUAD4 or SHELL4)
+ block_bc_flag.append(4)
+ block_bc.append(block)
+ bc[block]=4 #face has connectivity = 4
+ if name == self.topo: topography_face=block
+ else:
+ # block elements differ from HEX8/QUAD4/SHELL4
+ print '****************************************'
+ print 'block not properly defined:'
+ print ' name:',name
+ print ' type:',type
+ print
+ print 'please check your block definitions!'
+ print
+ print 'only supported types are:'
+ print ' HEX8 for volumes'
+ print ' QUAD4 for surface'
+ print ' SHELL4 for surface'
+ print '****************************************'
+ continue
+
+ nsets=cubit.get_nodeset_id_list()
+ if len(nsets) == 0: self.receivers=None
+ for nset in nsets:
+ name=cubit.get_exodus_entity_name('nodeset',nset)
+ if name == self.rec:
+ self.receivers=nset
+ else:
+ print 'nodeset '+name+' not defined'
+ self.receivers=None
+ try:
+ self.block_mat=block_mat
+ self.block_flag=block_flag
+ self.block_bc=block_bc
+ self.block_bc_flag=block_bc_flag
+ self.material=material
+ self.bc=bc
+ self.topography=topography_face
+ except:
+ print '****************************************'
+ print 'sorry, no blocks or blocks not properly defined'
+ print block_mat
+ print block_flag
+ print block_bc
+ print block_bc_flag
+ print material
+ print bc
+ print topography
+ print '****************************************'
+ def mat_parameter(self,properties):
+ #note: material property acoustic/elastic/poroelastic are defined by the block's name
+ print "#material properties:"
+ print properties
+ imaterial=properties[0]
+ flag=properties[1]
+ if flag > 0:
+ vel=properties[2]
+ if properties[3] is None and type(vel) != str:
+ # velocity model scales with given vp value
+ if vel >= 30:
+ m2km=1000.
+ else:
+ m2km=1.
+ vp=vel/m2km
+ rho=(1.6612*vp-0.472*vp**2+0.0671*vp**3-0.0043*vp**4+0.000106*vp**4)*m2km
+ txt='%1i %3i %20f %20f %20f %1i %1i\n' % (properties[0],properties[1],rho,vel,vel/(3**.5),0,0)
+ elif type(vel) != str:
+ # velocity model given as vp,vs,rho,..
+ #format nummaterials file: #material_domain_id #material_id #rho #vp #vs #Q_flag #anisotropy_flag
+ txt='%1i %3i %20f %20f %20f %2i %2i\n' % (properties[0],properties[1],properties[4], \
+ properties[2],properties[3],properties[5],properties[6])
+ else:
+ txt='%1i %3i %s \n' % (properties[0],properties[1],properties[2])
+ elif flag < 0:
+ if properties[2] == 'tomography':
+ txt='%1i %3i %s %s\n' % (properties[0],properties[1],properties[2],properties[3])
+ elif properties[2] == 'interface':
+ txt='%1i %3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],properties[3],\
+ properties[4],properties[5])
+ return txt
+ def nummaterial_write(self,nummaterial_name):
+ print 'Writing '+nummaterial_name+'.....'
+ nummaterial=open(nummaterial_name,'w')
+ for block in self.block_mat:
+ #name=cubit.get_exodus_entity_name('block',block)
+ nummaterial.write(self.mat_parameter(self.material[block]))
+ nummaterial.close()
+ print 'Ok'
+ def mesh_write(self,mesh_name):
+ meshfile=open(mesh_name,'w')
+ print 'Writing '+mesh_name+'.....'
+ num_elems=cubit.get_hex_count()
+ print ' number of elements:',str(num_elems)
+ meshfile.write(str(num_elems)+'\n')
+ num_write=0
+ for block,flag in zip(self.block_mat,self.block_flag):
+ #print block,flag
+ hexes=cubit.get_block_hexes(block)
+ #print len(hexes)
+ for hexa in hexes:
+ #print hexa
+ nodes=cubit.get_connectivity('Hex',hexa)
+ #nodes=self.jac_check(nodes) #is it valid for 3D? TODO
+ txt=('%10i ')% hexa
+ txt=txt+('%10i %10i %10i %10i %10i %10i %10i %10i\n')% nodes[:]
+ meshfile.write(txt)
+ meshfile.close()
+ print 'Ok'
+ def material_write(self,mat_name):
+ mat=open(mat_name,'w')
+ print 'Writing '+mat_name+'.....'
+ for block,flag in zip(self.block_mat,self.block_flag):
+ hexes=cubit.get_block_hexes(block)
+ for hexa in hexes:
+ mat.write(('%10i %10i\n') % (hexa,flag))
+ mat.close()
+ print 'Ok'
+ def nodescoord_write(self,nodecoord_name):
+ nodecoord=open(nodecoord_name,'w')
+ print 'Writing '+nodecoord_name+'.....'
+ node_list=cubit.parse_cubit_list('node','all')
+ num_nodes=len(node_list)
+ print ' number of nodes:',str(num_nodes)
+ nodecoord.write('%10i\n' % num_nodes)
+ #
+ for node in node_list:
+ x,y,z=cubit.get_nodal_coordinates(node)
+ txt=('%10i %20f %20f %20f\n') % (node,x,y,z)
+ nodecoord.write(txt)
+ nodecoord.close()
+ print 'Ok'
+ def free_write(self,freename=None):
+ # free surface
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ from sets import Set
+ normal=(0,0,1)
+ if not freename: freename=self.freename
+ # writes free surface file
+ print 'Writing '+freename+'.....'
+ freehex=open(freename,'w')
+ # searches block definition with name face_topo
+ for block,flag in zip(self.block_bc,self.block_bc_flag):
+ if block == self.topography:
+ name=cubit.get_exodus_entity_name('block',block)
+ print ' block name:',name,'id:',block
+ quads_all=cubit.get_block_faces(block)
+ print ' face = ',len(quads_all)
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ freehex.write('%10i\n' % len(quads_all))
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ #print f
+ nodes=cubit.get_connectivity('Face',f)
+ nodes_ok=self.normal_check(nodes,normal)
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ freehex.write(txt)
+ freehex.close()
+ print 'Ok'
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ def abs_write(self,absname=None):
+ # absorbing boundaries
+ import re
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ from sets import Set
+ if not absname: absname=self.absname
+ #
+ # loops through all block definitions
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for block,flag in zip(self.block_bc,self.block_bc_flag):
+ if block != self.topography:
+ name=cubit.get_exodus_entity_name('block',block)
+ print name,block
+ absflag=False
+ if re.search('xmin',name):
+ filename=absname+'_xmin'
+ normal=(-1,0,0)
+ elif re.search('xmax',name):
+ filename=absname+'_xmax'
+ normal=(1,0,0)
+ elif re.search('ymin',name):
+ filename=absname+'_ymin'
+ normal=(0,-1,0)
+ elif re.search('ymax',name):
+ filename=absname+'_ymax'
+ normal=(0,1,0)
+ elif re.search('bottom',name):
+ filename=absname+'_bottom'
+ normal=(0,0,-1)
+ elif re.search('abs',name):
+ print " ...face_abs - not used so far..."
+ continue
+ else:
+ continue
+ # opens file
+ print 'Writing '+filename+'.....'
+ abshex_local=open(filename,'w')
+ # gets face elements
+ quads_all=cubit.get_block_faces(block)
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ abshex_local.write('%10i\n' % len(quads_all))
+ #command = "group 'list_hex' add hex in face "+str(quads_all)
+ #command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
+ #cubit.cmd(command)
+ #group=cubit.get_id_from_name("list_hex")
+ #list_hex=cubit.get_group_hexes(group)
+ #command = "delete group "+ str(group)
+ #cubit.cmd(command)
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ nodes=cubit.get_connectivity('Face',f)
+ if not absflag:
+ # checks with specified normal
+ nodes_ok=self.normal_check(nodes,normal)
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ else:
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
+ abshex_local.write(txt)
+ # closes file
+ abshex_local.close()
+ print 'Ok'
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ def surface_write(self,pathdir=None):
+ # optional surfaces, e.g. moho_surface
+ # should be created like e.g.:
+ # > block 10 face in surface 2
+ # > block 10 name 'moho_surface'
+ import re
+ from sets import Set
+ for block in self.block_bc :
+ if block != self.topography:
+ name=cubit.get_exodus_entity_name('block',block)
+ # skips block names like face_abs**, face_topo**
+ if re.search('abs',name):
+ continue
+ elif re.search('topo',name):
+ continue
+ elif re.search('surface',name):
+ filename=pathdir+name+'_file'
+ else:
+ continue
+ # gets face elements
+ print ' surface block name: ',name,'id: ',block
+ quads_all=cubit.get_block_faces(block)
+ print ' face = ',len(quads_all)
+ if len(quads_all) == 0 :
+ continue
+ # writes out surface infos to file
+ print 'Writing '+filename+'.....'
+ surfhex_local=open(filename,'w')
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ # writes number of surface elements
+ surfhex_local.write('%10i\n' % len(quads_all))
+ # writes out element node ids
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ nodes=cubit.get_connectivity('Face',f)
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
+ surfhex_local.write(txt)
+ # closes file
+ surfhex_local.close()
+ print 'Ok'
+ def rec_write(self,recname):
+ print 'Writing '+self.recname+'.....'
+ recfile=open(self.recname,'w')
+ nodes=cubit.get_nodeset_nodes(self.receivers)
+ for i,n in enumerate(nodes):
+ x,y,z=cubit.get_nodal_coordinates(n)
+ recfile.write('ST%i XX %20f %20f 0.0 0.0 \n' % (i,x,z))
+ recfile.close()
+ print 'Ok'
+ def write(self,path=''):
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ if len(path) != 0:
+ if path[-1] != '/': path=path+'/'
+ # mesh file
+ self.mesh_write(path+self.mesh_name)
+ # mesh material
+ self.material_write(path+self.material_name)
+ # mesh coordinates
+ self.nodescoord_write(path+self.nodecoord_name)
+ # material definitions
+ self.nummaterial_write(path+self.nummaterial_name)
+ # free surface: face_top
+ self.free_write(path+self.freename)
+ # absorbing surfaces: abs_***
+ self.abs_write(path+self.absname)
+ # any other surfaces: ***surface***
+ self.surface_write(path)
+ # receivers
+ if self.receivers: self.rec_write(path+self.recname)
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+
+def export2SESAME(path_exporting_mesh_SPECFEM3D_SESAME):
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ sem_mesh=mesh()
+ sem_mesh.write(path=path_exporting_mesh_SPECFEM3D_SESAME)
+
+
+if __name__ == '__main__':
+ path='MESH/'
+ export2SESAME(path)
+
+# call by:
+# import cubit2specfem3d
+# cubit2specfem3d.export2SESAME('MESH')
Added: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/run_boundary_definition.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/run_boundary_definition.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/run_boundary_definition.py 2010-08-16 20:37:29 UTC (rev 17093)
@@ -0,0 +1,18 @@
+#!python
+#!/usr/bin/env python
+
+import cubit
+import boundary_definition
+import cubit2specfem3d
+
+import os
+import sys
+
+
+###### This is boundary_definition.py of GEOCUBIT
+#..... which extracts the bounding faces and defines them into blocks
+reload(boundary_definition)
+boundary_definition.entities=['face']
+boundary_definition.define_bc(boundary_definition.entities,parallel=True)
+
+
Property changes on: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/run_boundary_definition.py
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/run_cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/run_cubit2specfem3d.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/run_cubit2specfem3d.py 2010-08-16 20:37:29 UTC (rev 17093)
@@ -0,0 +1,26 @@
+#!python
+#!/usr/bin/env python
+
+import cubit
+import boundary_definition
+import cubit2specfem3d
+
+import os
+import sys
+
+
+###### This is boundary_definition.py of GEOCUBIT
+#..... which extracts the bounding faces and defines them into blocks
+#reload(boundary_definition)
+#boundary_definition.entities=['face']
+#boundary_definition.define_bc(boundary_definition.entities,parallel=True)
+
+
+#### Export to SESAME format using cubit2specfem3d.py of GEOCUBIT
+os.system('mkdir -p MESH')
+
+reload(cubit2specfem3d)
+cubit2specfem3d.export2SESAME('MESH')
+
+# all files needed by SCOTCH are now in directory MESH
+
Property changes on: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/run_cubit2specfem3d.py
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/tomoblock_mesh.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/tomoblock_mesh.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/tomoblock_mesh.py 2010-08-16 20:37:29 UTC (rev 17093)
@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+
+import cubit
+import boundary_definition
+import cubit2specfem3d
+
+import os
+import sys
+
+cubit.cmd('reset')
+cubit.cmd('brick x 134000 y 134000 z 60000')
+cubit.cmd('volume 1 move x 67000 y 67000 z -30000')
+
+
+# Meshing the volumes
+elementsize = 3750.0
+
+cubit.cmd('volume 1 size '+str(elementsize))
+cubit.cmd('mesh volume 1')
+
+
+#### End of meshing
+
+###### This is boundary_definition.py of GEOCUBIT
+#..... which extracts the bounding faces and defines them into blocks
+boundary_definition.entities=['face']
+boundary_definition.define_bc(boundary_definition.entities,parallel=True)
+
+#### Define material properties for the 3 volumes ################
+cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
+
+
+cubit.cmd('block 1 name "elastic tomography_model.xyz 1" ') # elastic material region
+cubit.cmd('block 1 attribute count 2')
+cubit.cmd('block 1 attribute index 1 -1') # flag for material: -1 for 1. undefined material
+cubit.cmd('block 1 attribute index 2 2') # flag for tomographic model
+
+
+#cubit.cmd('block 1 name "acoustic tomographic 1" ') # acoustic material region
+#cubit.cmd('block 1 attribute count 2')
+#cubit.cmd('block 1 attribute index 1 -1') # material 1
+#cubit.cmd('block 1 attribute index 2 2') # tomographic model flag
+
+
+cubit.cmd('export mesh "top.e" dimension 3 overwrite')
+cubit.cmd('save as "meshing.cub" overwrite')
+
+#### Export to SESAME format using cubit2specfem3d.py of GEOCUBIT
+
+os.system('mkdir -p MESH')
+cubit2specfem3d.export2SESAME('MESH')
+
+# all files needed by SCOTCH are now in directory MESH
Property changes on: seismo/3D/SPECFEM3D/trunk/EXAMPLES/tomographic_model/tomoblock_mesh.py
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/SPECFEM3D/trunk/EXAMPLES/waterlayered_halfspace/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/EXAMPLES/waterlayered_halfspace/cubit2specfem3d.py 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/EXAMPLES/waterlayered_halfspace/cubit2specfem3d.py 2010-08-16 20:37:29 UTC (rev 17093)
@@ -539,7 +539,7 @@
meshfile=open(mesh_name,'w')
print 'Writing '+mesh_name+'.....'
num_elems=cubit.get_hex_count()
- print num_elems
+ print ' number of elements:',str(num_elems)
meshfile.write(str(num_elems)+'\n')
num_write=0
for block,flag in zip(self.block_mat,self.block_flag):
@@ -554,7 +554,7 @@
txt=txt+('%10i %10i %10i %10i %10i %10i %10i %10i\n')% nodes[:]
meshfile.write(txt)
meshfile.close()
- print 'Ok',str(num_elems)
+ print 'Ok'
def material_write(self,mat_name):
mat=open(mat_name,'w')
print 'Writing '+mat_name+'.....'
@@ -569,6 +569,7 @@
print 'Writing '+nodecoord_name+'.....'
node_list=cubit.parse_cubit_list('node','all')
num_nodes=len(node_list)
+ print ' number of nodes:',str(num_nodes)
nodecoord.write('%10i\n' % num_nodes)
#
for node in node_list:
Modified: seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -114,18 +114,18 @@
f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
t0 = 1.2d0/f0
- if (it == 1 .and. myrank == 0) then
- print *,'using a source of dominant frequency ',f0
- print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
- print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
- endif
+ !if (it == 1 .and. myrank == 0) then
+ ! write(IMAIN,*) 'using a source of dominant frequency ',f0
+ ! write(IMAIN,*) 'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+ ! write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+ !endif
! gaussian source time function
!stf_used = comp_source_time_function(dble(it-1)*DT-t0-t_cmt(isource),hdur_gaussian(isource))
! we use nu_source(:,3) here because we want a source normal to the surface.
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = 1.d10 * comp_source_time_function_rickr(dble(it-1)*DT-t0-t_cmt(isource),f0)
+ stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(it-1)*DT-t0-t_cmt(isource),f0)
! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
! the sign is negative because pressure p = - Chi_dot_dot therefore we need
@@ -256,18 +256,18 @@
f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
t0 = 1.2d0/f0
- if (it == 1 .and. myrank == 0) then
- print *,'using a source of dominant frequency ',f0
- print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
- print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
- endif
+ !if (it == 1 .and. myrank == 0) then
+ ! write(IMAIN,*) 'using a source of dominant frequency ',f0
+ ! write(IMAIN,*) 'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+ ! write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+ !endif
! gaussian source time function
!stf_used = comp_source_time_function(dble(it-1)*DT-t0-t_cmt(isource),hdur_gaussian(isource))
! we use nu_source(:,3) here because we want a source normal to the surface.
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = 1.d10 * comp_source_time_function_rickr(dble(NSTEP-it-1)*DT-t0-t_cmt(isource),f0)
+ stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(NSTEP-it-1)*DT-t0-t_cmt(isource),f0)
! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
! the sign is negative because pressure p = - Chi_dot_dot therefore we need
Modified: seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -110,16 +110,16 @@
f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
t0 = 1.2d0/f0
- if (it == 1 .and. myrank == 0) then
- print *,'using a source of dominant frequency ',f0
- print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
- print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
- endif
+ !if (it == 1 .and. myrank == 0) then
+ ! write(IMAIN,*) 'using a source of dominant frequency ',f0
+ ! write(IMAIN,*) 'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+ ! write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+ !endif
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = 1.d10 * comp_source_time_function_rickr(dble(it-1)*DT-t0-t_cmt(isource),f0)
+ stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(it-1)*DT-t0-t_cmt(isource),f0)
- ! we use nu_source(:,3) here because we want a source normal to the surface.
+ ! we use nu_source(:,3) here because we want a source normal to the surface (z-direction).
accel(:,iglob) = accel(:,iglob) &
+ sngl( nu_source(:,3,isource) ) * stf_used
@@ -231,14 +231,14 @@
f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
t0 = 1.2d0/f0
- if (it == 1 .and. myrank == 0) then
- print *,'using a source of dominant frequency ',f0
- print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
- print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
- endif
+ !if (it == 1 .and. myrank == 0) then
+ ! write(IMAIN,*) 'using a source of dominant frequency ',f0
+ ! write(IMAIN,*) 'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+ ! write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+ !endif
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = 1.d10 * comp_source_time_function_rickr(dble(NSTEP-it-1)*DT-t0-t_cmt(isource),f0)
+ stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(NSTEP-it-1)*DT-t0-t_cmt(isource),f0)
! we use nu_source(:,3) here because we want a source normal to the surface.
! note: time step is now at NSTEP-it
Modified: seismo/3D/SPECFEM3D/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/create_regions_mesh.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/create_regions_mesh.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -361,7 +361,8 @@
ispec_is_acoustic,ispec_is_elastic,ispec_is_poroelastic)
! computes the approximate amount of static memory needed to run the solver
- call memory_eval(nspec,nglob,maxval(nibool_interfaces_ext_mesh),num_interfaces_ext_mesh,static_memory_size)
+ call memory_eval(nspec,nglob,maxval(nibool_interfaces_ext_mesh),num_interfaces_ext_mesh, &
+ OCEANS,static_memory_size)
call max_all_dp(static_memory_size, max_static_memory_size)
! checks the mesh, stability and resolved period
Modified: seismo/3D/SPECFEM3D/trunk/decompose_mesh_SCOTCH/Makefile
===================================================================
--- seismo/3D/SPECFEM3D/trunk/decompose_mesh_SCOTCH/Makefile 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/decompose_mesh_SCOTCH/Makefile 2010-08-16 20:37:29 UTC (rev 17093)
@@ -3,7 +3,7 @@
#############################################################
## modify to match your compiler defaults
## (which were used to compile SCOTCH libraries from below as well)
-F90 = ifort # use -warn
+F90 = ifort # use -g -traceback -check bounds -warn
#F90 = gfortran # use -Wall
## modify to match your library paths
Modified: seismo/3D/SPECFEM3D/trunk/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -88,7 +88,9 @@
! reads in mesh files
!----------------------------------------------------------------------------------------------
subroutine read_mesh_files
-
+ implicit none
+ character(len=256) :: line
+
! sets number of nodes per element
ngnod = esize
@@ -192,10 +194,12 @@
print *,'materials:'
! counts materials (defined/undefined)
do while (ierr == 0)
- print*, ' num_mat = ',num_mat
- if(num_mat /= -1) then
+ print*, ' num_mat = ',num_mat
+ if(num_mat > 0 ) then
+ ! positive materials_id: velocity values will be defined
count_def_mat = count_def_mat + 1
else
+ ! negative materials_id: undefined material properties yet
count_undef_mat = count_undef_mat + 1
end if
read(98,*,iostat=ierr) idummy,num_mat
@@ -214,6 +218,10 @@
! reads in defined material properties
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/nummaterial_velocity_file', &
status='old', form='formatted')
+
+ ! note: entries in nummaterial_velocity_file must be sorted to list all
+ ! defined materials (material_id > 0) first, and afterwards list all
+ ! undefined materials (material_id < 0 )
do imat=1,count_def_mat
! material definitions
!
@@ -221,6 +229,9 @@
! #(6) material_domain_id #(0) material_id #(1) rho #(2) vp #(3) vs #(4) Q_flag #(5) anisotropy_flag
!
read(98,*) idomain_id,num_mat,rho,vp,vs,q_flag,aniso_flag
+
+ if(num_mat < 1 .or. num_mat > count_def_mat) stop "ERROR : Invalid nummaterial_velocity_file file."
+
!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)
mat_prop(1,num_mat) = rho
@@ -228,9 +239,7 @@
mat_prop(3,num_mat) = vs
mat_prop(4,num_mat) = q_flag
mat_prop(5,num_mat) = aniso_flag
- mat_prop(6,num_mat) = idomain_id
-
- if(num_mat < 0 .or. num_mat > count_def_mat) stop "ERROR : Invalid nummaterial_velocity_file file."
+ mat_prop(6,num_mat) = idomain_id
!checks attenuation flag with integer range as defined in constants.h like IATTENUATION_SEDIMENTS_40, ....
if( int(mat_prop(4,num_mat)) > 13 ) then
@@ -239,11 +248,97 @@
end do
! reads in undefined material properties
do imat=1,count_undef_mat
- read(98,'(6A30)') undef_mat_prop(6,imat),undef_mat_prop(1,imat),undef_mat_prop(2,imat),&
+ ! undefined materials: have to be listed in decreasing order of material_id (start with -1, -2, etc...)
+ ! format:
+ ! - for interfaces
+ ! #material_domain_id #material_id(<0) #type_name (="interface") #material_id_for_material_below #material_id_for_material_above
+ ! example: 2 -1 interface 1 2
+ ! - for tomography models
+ ! #material_domain_id #material_id (<0) #type_name (="tomography") #block_name
+ ! example: 2 -1 tomography elastic tomography_model.xyz 1
+ read(98,'(A256)') line
+
+ ! checks if interface or tomography definition
+ read(line,*) undef_mat_prop(6,imat),undef_mat_prop(1,imat),undef_mat_prop(2,imat)
+ if( trim(undef_mat_prop(2,imat)) == 'interface' ) then
+ ! line will have 5 arguments, e.g.: 2 -1 interface 1 2
+ read(line,*) undef_mat_prop(6,imat),undef_mat_prop(1,imat),undef_mat_prop(2,imat),&
+ undef_mat_prop(3,imat),undef_mat_prop(4,imat)
+ undef_mat_prop(5,imat) = "1" ! dummy value
+ else if( trim(undef_mat_prop(2,imat)) == 'tomography' ) then
+ ! line will have 6 arguments, e.g.: 2 -1 tomography elastic tomography_model.xyz 1
+ read(line,*) undef_mat_prop(6,imat),undef_mat_prop(1,imat),undef_mat_prop(2,imat),&
undef_mat_prop(3,imat),undef_mat_prop(4,imat),undef_mat_prop(5,imat)
+ else
+ stop "ERROR: invalid line in nummaterial_velocity_file for undefined material"
+ endif
+
+ !read(98,*) undef_mat_prop(6,imat),undef_mat_prop(1,imat),undef_mat_prop(2,imat),&
+ ! undef_mat_prop(3,imat),undef_mat_prop(4,imat),undef_mat_prop(5,imat)
+
+ ! debug output
+ !print*,'properties:'
+ !print*,undef_mat_prop(:,imat)
+ !print*
+
+ ! checks material_id
+ read(undef_mat_prop(1,imat),*) num_mat
+ !print *,'material_id: ',num_mat
+ if(num_mat > 0 .or. -num_mat > count_undef_mat) stop "ERROR : Invalid nummaterial_velocity_file for undefined materials."
+ if(num_mat /= -imat) stop "ERROR : Invalid material_id in nummaterial_velocity_file for undefined materials."
+
+ ! checks interface: flag_down/flag_up
+ if( trim(undef_mat_prop(2,imat)) == 'interface' ) then
+ ! flag_down
+ read( undef_mat_prop(3,imat),*) num_mat
+ if( num_mat > 0 ) then
+ ! must point to a defined material
+ if( num_mat > count_def_mat) stop "ERROR: invalid flag_down in interface definition in nummaterial_velocity_file"
+ else
+ ! must point to an undefined material
+ if( -num_mat > count_undef_mat) stop "ERROR: invalid flag_down in interface definition in nummaterial_velocity_file"
+ endif
+ ! flag_up
+ read( undef_mat_prop(4,imat),*) num_mat
+ if( num_mat > 0 ) then
+ ! must point to a defined material
+ if( num_mat > count_def_mat) stop "ERROR: invalid flag_up in interface definition in nummaterial_velocity_file"
+ else
+ ! must point to an undefined material
+ if( -num_mat > count_undef_mat) stop "ERROR: invalid flag_up in interface definition in nummaterial_velocity_file"
+ endif
+ endif
end do
close(98)
+
+ ! TODO:
+ ! must be changed, if mat(1,i) < 0 1 == interface , 2 == tomography
+ do ispec=1,nspec
+ ! get material_id
+ num_mat = mat(1,ispec)
+ if( num_mat < 0 ) then
+ ! finds undefined material property
+ do imat=1,count_undef_mat
+ if( -imat == num_mat ) then
+ ! interface
+ if( trim(undef_mat_prop(2,imat)) == 'interface' ) then
+ mat(2,ispec) = 1
+ ! tomography
+ elseif( trim(undef_mat_prop(2,imat)) == 'tomography' ) then
+ mat(2,ispec) = 2
+ else
+ ! shouldn't encounter this case
+ stop "error undefined material: type name not recognized"
+ endif
+ endif
+ enddo
+ else
+ ! ispec belongs to a defined material
+ mat(2,ispec) = 0
+ endif
+ enddo
+
! reads in absorbing boundary files
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_xmin', &
status='old', form='formatted',iostat=ierr)
@@ -452,19 +547,44 @@
elmnts_load(:) = 1
! in case of acoustic/elastic simulation, weights elements accordingly
- call acoustic_elastic_load(elmnts_load,nspec,count_def_mat,mat(1,:),mat_prop)
+ call acoustic_elastic_load(elmnts_load,nspec,count_def_mat,count_undef_mat, &
+ mat(1,:),mat_prop,undef_mat_prop)
! SCOTCH partitioning
+
+ ! we use default strategy for partitioning, thus omit specifing explicit strategy .
+
+ ! workflow preferred by F. Pellegrini (SCOTCH):
+ !!This comes from the fact that, in version 5.1.8, the name
+ !!for the "recursive bisection" method has changed from "b"
+ !!("bipartitioning") to "r" ("recursive").
+ !!
+ !!As a general rule, do not try to set up strategies by
+ !!yourself. The default strategy in Scotch will most probably
+ !!provide better results. To use it, just call:
+ !!
+ !!SCOTCHFstratInit (),
+ !!
+ !!and use this "empty" strategy in the mapping routine
+ !!(consequently, no call to SCOTCHFstratGraphMap () is
+ !!required).
+ !!
+ !!This will make you independent from further changes (most
+ !!probably improvements !;-) ) in the strategy syntax.
+ !!And you should see an improvement in performance, too,
+ !!as your hand-made strategy did not make use of the
+ !!multi-level framework.
+
call scotchfstratinit (scotchstrat(1), ierr)
if (ierr /= 0) then
stop 'ERROR : MAIN : Cannot initialize strat'
- endif
+ endif
+
+ !call scotchfstratgraphmap (scotchstrat(1), trim(scotch_strategy), ierr)
+ ! if (ierr /= 0) then
+ ! stop 'ERROR : MAIN : Cannot build strat'
+ !endif
- call scotchfstratgraphmap (scotchstrat(1), trim(scotch_strategy), ierr)
- if (ierr /= 0) then
- stop 'ERROR : MAIN : Cannot build strat'
- endif
-
call scotchfgraphinit (scotchgraph (1), ierr)
if (ierr /= 0) then
stop 'ERROR : MAIN : Cannot initialize graph'
Modified: seismo/3D/SPECFEM3D/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -581,9 +581,9 @@
mat_prop(4,i), mat_prop(5,i), mat_prop(6,i)
end do
do i = 1, count_undef_mat
- write(IIN_database,*) trim(undef_mat_prop(1,i)),trim(undef_mat_prop(2,i)), &
- trim(undef_mat_prop(3,i)),trim(undef_mat_prop(4,i)), &
- trim(undef_mat_prop(5,i)),trim(undef_mat_prop(6,i))
+ write(IIN_database,*) trim(undef_mat_prop(1,i)),' ',trim(undef_mat_prop(2,i)),' ', &
+ trim(undef_mat_prop(3,i)),' ',trim(undef_mat_prop(4,i)),' ', &
+ trim(undef_mat_prop(5,i)),' ',trim(undef_mat_prop(6,i))
end do
end subroutine write_material_properties_database
@@ -1159,7 +1159,8 @@
! expensive calculations in specfem simulations
!--------------------------------------------------
- subroutine acoustic_elastic_load (elmnts_load,nelmnts,nb_materials,num_material,mat_prop)
+ subroutine acoustic_elastic_load (elmnts_load,nelmnts,count_def_mat,count_undef_mat, &
+ num_material,mat_prop,undef_mat_prop)
!
! note:
! acoustic material = domainID 1 (stored in mat_prop(6,..) )
@@ -1168,33 +1169,51 @@
implicit none
integer(long),intent(in) :: nelmnts
- integer, intent(in) :: nb_materials
+ integer, intent(in) :: count_def_mat,count_undef_mat
! load weights
integer,dimension(1:nelmnts),intent(out) :: elmnts_load
! materials
integer, dimension(1:nelmnts), intent(in) :: num_material
- double precision, dimension(6,nb_materials),intent(in) :: mat_prop
+ double precision, dimension(6,count_def_mat),intent(in) :: mat_prop
+ character (len=30), dimension(6,count_undef_mat),intent(in) :: undef_mat_prop
! local parameters
- logical, dimension(nb_materials) :: is_acoustic, is_elastic
- integer :: i,el
-
- ! sets acoustic/elastic flags for materials
+ logical, dimension(-count_undef_mat:count_def_mat) :: is_acoustic, is_elastic
+ integer :: i,el,idomain_id
+
+ ! initializes flags
is_acoustic(:) = .false.
is_elastic(:) = .false.
- do i = 1, nb_materials
- ! acoustic material has idomain_id 1
- if (mat_prop(6,i) == 1 ) then
+
+ ! sets acoustic/elastic flags for defined materials
+ do i = 1, count_def_mat
+ idomain_id = mat_prop(6,i)
+ ! acoustic material has idomain_id 1
+ if (idomain_id == 1 ) then
is_acoustic(i) = .true.
endif
! elastic material has idomain_id 2
- if (mat_prop(6,i) == 2 ) then
+ if (idomain_id == 2 ) then
is_elastic(i) = .true.
endif
enddo
+ ! sets acoustic/elastic flags for undefined materials
+ do i = 1, count_undef_mat
+ read(undef_mat_prop(6,i),'(i)') idomain_id
+ ! acoustic material has idomain_id 1
+ if (idomain_id == 1 ) then
+ is_acoustic(-i) = .true.
+ endif
+ ! elastic material has idomain_id 2
+ if (idomain_id == 2 ) then
+ is_elastic(-i) = .true.
+ endif
+ enddo
+
+
! sets weights for elements
do el = 0, nelmnts-1
! acoustic element (cheap)
Modified: seismo/3D/SPECFEM3D/trunk/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/generate_databases.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/generate_databases.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -373,8 +373,10 @@
! check that the code is running with the requested nb of processes
if(sizeprocs /= NPROC) then
- write(IMAIN,*) 'error: number of processors supposed to run on: ',NPROC
- write(IMAIN,*) 'error: number of processors actually run on: ',sizeprocs
+ if( myrank == 0 ) then
+ write(IMAIN,*) 'error: number of processors supposed to run on: ',NPROC
+ write(IMAIN,*) 'error: number of processors actually run on: ',sizeprocs
+ endif
call exit_MPI(myrank,'wrong number of MPI processes')
endif
@@ -382,7 +384,7 @@
! just to be sure for now..
if( ABSORBING_CONDITIONS ) then
if( NGLLX /= NGLLY .and. NGLLY /= NGLLZ ) &
- stop 'must have NGLLX = NGLLY = NGLLZ for external meshes'
+ call exit_MPI(myrank,'must have NGLLX = NGLLY = NGLLZ for external meshes')
endif
! info about external mesh simulation
@@ -555,7 +557,6 @@
! read materials' physical properties
read(IIN,*) nmat_ext_mesh, nundefMat_ext_mesh
allocate(materials_ext_mesh(6,nmat_ext_mesh))
- allocate(undef_mat_prop(6,nundefMat_ext_mesh))
do imat = 1, nmat_ext_mesh
! format: #(1) rho #(2) vp #(3) vs #(4) Q_flag #(5) anisotropy_flag #(6) material_domain_id
read(IIN,*) materials_ext_mesh(1,imat), materials_ext_mesh(2,imat), materials_ext_mesh(3,imat), &
@@ -571,9 +572,18 @@
endif
call sync_all()
+ allocate(undef_mat_prop(6,nundefMat_ext_mesh))
do imat = 1, nundefMat_ext_mesh
+ ! format example tomography:
+ ! -1 tomography elastic tomography_model.xyz 1 2
+ ! format example interface:
+ ! -1 interface 14 15 1 2
read(IIN,*) undef_mat_prop(1,imat),undef_mat_prop(2,imat),undef_mat_prop(3,imat),undef_mat_prop(4,imat), &
undef_mat_prop(5,imat), undef_mat_prop(6,imat)
+
+ ! output debug
+ !print*,'undefined materials:'
+ !print*,undef_mat_prop(:,imat)
end do
if(myrank == 0) then
@@ -593,6 +603,10 @@
read(IIN,*) dummy_elmnt, mat_ext_mesh(1,ispec),mat_ext_mesh(2,ispec), &
elmnts_ext_mesh(1,ispec), elmnts_ext_mesh(2,ispec), elmnts_ext_mesh(3,ispec), elmnts_ext_mesh(4,ispec), &
elmnts_ext_mesh(5,ispec), elmnts_ext_mesh(6,ispec), elmnts_ext_mesh(7,ispec), elmnts_ext_mesh(8,ispec)
+
+ ! check debug
+ if( dummy_elmnt /= ispec) stop "error ispec order in materials file"
+
enddo
NSPEC_AB = nelmnts_ext_mesh
Modified: seismo/3D/SPECFEM3D/trunk/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/get_model.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/get_model.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -123,15 +123,16 @@
do iundef = 1,nundefMat_ext_mesh
if(trim(undef_mat_prop(2,iundef)) == 'interface') then
- read(undef_mat_prop(4,iundef),'(1i3)') flag_below
- read(undef_mat_prop(5,iundef),'(1i3)') flag_above
+ read(undef_mat_prop(3,iundef),'(1i3)') flag_below
+ read(undef_mat_prop(4,iundef),'(1i3)') flag_above
endif
enddo
! see file model_interface_bedrock.f90: routine interface()
!call interface(iflag,flag_below,flag_above,ispec,nspec,i,j,k,xstore,ystore,zstore,ibedrock)
-
- iflag = 1
+
+ ! dummy: takes 1. defined material
+ iflag = 1
rho = materials_ext_mesh(1,iflag)
vp = materials_ext_mesh(2,iflag)
vs = materials_ext_mesh(3,iflag)
@@ -145,7 +146,7 @@
iflag_aniso = materials_ext_mesh(5,iflag)
idomain_id = materials_ext_mesh(6,iflag)
- else if ( imaterial_id < 0 ) then
+ else if ( mat_ext_mesh(2,ispec) == 2 ) then
! material definition undefined, uses definition from tomography model
! GLL point location
@@ -160,7 +161,12 @@
iflag_atten = 1 ! attenuation: would use IATTENUATION_SEDIMENTS_40
iflag_aniso = 0 ! no anisotropy
- idomain_id = 2 ! elastic domain
+
+ ! sets acoustic/elastic domain as given in materials properties
+ iundef = - imaterial_id ! iundef must be positive
+ read(undef_mat_prop(6,iundef),*) idomain_id
+ ! or
+ !idomain_id = IDOMAIN_ELASTIC ! forces to be elastic domain
else
Modified: seismo/3D/SPECFEM3D/trunk/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/initialize_simulation.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/initialize_simulation.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -146,10 +146,7 @@
allocate(ispec_is_acoustic(NSPEC_AB))
allocate(ispec_is_elastic(NSPEC_AB))
allocate(ispec_is_poroelastic(NSPEC_AB))
-
- ! ocean mass matrix
- allocate(rmass_ocean_load(NGLOB_AB))
-
+
! initializes adjoint simulations
call initialize_simulation_adjoint()
Modified: seismo/3D/SPECFEM3D/trunk/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/locate_source.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/locate_source.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -54,7 +54,7 @@
integer myrank
-! arrays containing coordinates of the points
+ ! arrays containing coordinates of the points
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
logical, dimension(NSPEC_AB) :: ispec_is_acoustic,ispec_is_elastic
@@ -66,7 +66,8 @@
integer iprocloop
- integer i,j,k,ispec,iglob,iglob_selected,inode,iface,isource,imin,imax,jmin,jmax,kmin,kmax,igll,jgll,kgll
+ integer i,j,k,ispec,iglob,iglob_selected,inode,iface,isource
+ integer imin,imax,jmin,jmax,kmin,kmax,igll,jgll,kgll
integer iselected,jselected,iface_selected,iadjust,jadjust
integer iproc(1)
@@ -121,7 +122,8 @@
double precision, dimension(:), allocatable :: tmp_local
double precision, dimension(:,:),allocatable :: tmp_all_local
- double precision hdur(NSOURCES) !, hdur_gaussian(NSOURCES) !, t0
+ double precision hdur(NSOURCES)
+ double precision :: f0,t0_ricker
double precision, dimension(NSOURCES) :: Mxx,Myy,Mzz,Mxy,Mxz,Myz
double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
@@ -150,15 +152,9 @@
integer ix_initial_guess_source,iy_initial_guess_source,iz_initial_guess_source
- ! for calculation of source time function
- !integer it
- !double precision time_source
- !double precision, external :: comp_source_time_function
-
integer, dimension(NSOURCES) :: idomain
integer, dimension(NGATHER_SOURCES,0:NPROC-1) :: idomain_all
-
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
@@ -172,12 +168,6 @@
if(hdur(isource) < 5. * DT) hdur(isource) = 5. * DT
enddo
- ! convert the half duration for triangle STF to the one for gaussian STF
- !hdur_gaussian = hdur/SOURCE_DECAY_MIMIC_TRIANGLE
-
- ! define t0 as the earliest start time
- !t0 = - 1.5d0 * minval(t_cmt-hdur)
-
! define topology of the control element
call usual_hex_nodes(iaddx,iaddy,iaddz)
@@ -213,77 +203,77 @@
! set distance to huge initial value
distmin = HUGEVAL
if(num_free_surface_faces > 0) then
- iglob_selected = 1
- ! loop only on points inside the element
- ! exclude edges to ensure this point is not shared with other elements
- imin = 2
- imax = NGLLX - 1
+ iglob_selected = 1
+ ! loop only on points inside the element
+ ! exclude edges to ensure this point is not shared with other elements
+ imin = 2
+ imax = NGLLX - 1
- jmin = 2
- jmax = NGLLY - 1
- do iface=1,num_free_surface_faces
- do j=jmin,jmax
- do i=imin,imax
+ jmin = 2
+ jmax = NGLLY - 1
+ do iface=1,num_free_surface_faces
+ do j=jmin,jmax
+ do i=imin,imax
- ispec = free_surface_ispec(iface)
- igll = free_surface_ijk(1,(j-1)*NGLLY+i,iface)
- jgll = free_surface_ijk(2,(j-1)*NGLLY+i,iface)
- kgll = free_surface_ijk(3,(j-1)*NGLLY+i,iface)
- iglob = ibool(igll,jgll,kgll,ispec)
+ ispec = free_surface_ispec(iface)
+ igll = free_surface_ijk(1,(j-1)*NGLLY+i,iface)
+ jgll = free_surface_ijk(2,(j-1)*NGLLY+i,iface)
+ kgll = free_surface_ijk(3,(j-1)*NGLLY+i,iface)
+ iglob = ibool(igll,jgll,kgll,ispec)
- ! keep this point if it is closer to the receiver
- dist = dsqrt((utm_x_source(isource)-dble(xstore(iglob)))**2 + &
+ ! keep this point if it is closer to the receiver
+ dist = dsqrt((utm_x_source(isource)-dble(xstore(iglob)))**2 + &
(utm_y_source(isource)-dble(ystore(iglob)))**2)
- if(dist < distmin) then
- distmin = dist
- iglob_selected = iglob
- iface_selected = iface
- iselected = i
- jselected = j
- altitude_source(1) = zstore(iglob_selected)
- endif
- enddo
+ if(dist < distmin) then
+ distmin = dist
+ iglob_selected = iglob
+ iface_selected = iface
+ iselected = i
+ jselected = j
+ altitude_source(1) = zstore(iglob_selected)
+ endif
enddo
- ! end of loop on all the elements on the free surface
- end do
-! weighted mean at current point of topography elevation of the four closest nodes
-! set distance to huge initial value
- distmin = HUGEVAL
- do j=jselected,jselected+1
- do i=iselected,iselected+1
- inode = 1
- do jadjust=0,1
- do iadjust= 0,1
- ispec = free_surface_ispec(iface_selected)
- igll = free_surface_ijk(1,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
- jgll = free_surface_ijk(2,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
- kgll = free_surface_ijk(3,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
- iglob = ibool(igll,jgll,kgll,ispec)
+ enddo
+ ! end of loop on all the elements on the free surface
+ end do
+ ! weighted mean at current point of topography elevation of the four closest nodes
+ ! set distance to huge initial value
+ distmin = HUGEVAL
+ do j=jselected,jselected+1
+ do i=iselected,iselected+1
+ inode = 1
+ do jadjust=0,1
+ do iadjust= 0,1
+ ispec = free_surface_ispec(iface_selected)
+ igll = free_surface_ijk(1,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
+ jgll = free_surface_ijk(2,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
+ kgll = free_surface_ijk(3,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
+ iglob = ibool(igll,jgll,kgll,ispec)
- elevation_node(inode) = zstore(iglob)
- dist_node(inode) = dsqrt((utm_x_source(isource)-dble(xstore(iglob)))**2 + &
+ elevation_node(inode) = zstore(iglob)
+ dist_node(inode) = dsqrt((utm_x_source(isource)-dble(xstore(iglob)))**2 + &
(utm_y_source(isource)-dble(ystore(iglob)))**2)
- inode = inode + 1
- end do
- end do
- dist = sum(dist_node)
- if(dist < distmin) then
- distmin = dist
- altitude_source(1) = (dist_node(1)/dist)*elevation_node(1) + &
+ inode = inode + 1
+ end do
+ end do
+ dist = sum(dist_node)
+ if(dist < distmin) then
+ distmin = dist
+ altitude_source(1) = (dist_node(1)/dist)*elevation_node(1) + &
(dist_node(2)/dist)*elevation_node(2) + &
(dist_node(3)/dist)*elevation_node(3) + &
(dist_node(4)/dist)*elevation_node(4)
- endif
- end do
- end do
+ endif
+ end do
+ end do
end if
! MPI communications to determine the best slice
distmin_ele(1)= distmin
call gather_all_dp(distmin_ele,1,distmin_ele_all,1,NPROC)
call gather_all_dp(altitude_source,1,elevation_all,1,NPROC)
if(myrank == 0) then
- iproc = minloc(distmin_ele_all)
- altitude_source(1) = elevation_all(iproc(1))
+ iproc = minloc(distmin_ele_all)
+ altitude_source(1) = elevation_all(iproc(1))
end if
call bcast_all_dp(altitude_source,1)
elevation(isource) = altitude_source(1)
@@ -304,7 +294,8 @@
x_target_source = utm_x_source(isource)
y_target_source = utm_y_source(isource)
- !z_target_source = depth(isource)
+
+ ! depth in CMTSOLUTION given in km
z_target_source = - depth(isource)*1000.0d0 + elevation(isource)
! set distance to huge initial value
@@ -314,7 +305,6 @@
do ispec=1,NSPEC_AB
-
! define the interval in which we look for points
if(USE_FORCE_POINT_SOURCE) then
imin = 1
@@ -351,13 +341,13 @@
endif
endif
- ! keep this point if it is closer to the source
- dist=dsqrt((x_target_source-dble(xstore(iglob)))**2 &
+ ! keep this point if it is closer to the source
+ dist = dsqrt((x_target_source-dble(xstore(iglob)))**2 &
+(y_target_source-dble(ystore(iglob)))**2 &
+(z_target_source-dble(zstore(iglob)))**2)
if(dist < distmin) then
- distmin=dist
- ispec_selected_source(isource)=ispec
+ distmin = dist
+ ispec_selected_source(isource) = ispec
ix_initial_guess_source = i
iy_initial_guess_source = j
iz_initial_guess_source = k
@@ -389,9 +379,9 @@
! sets whether acoustic (1) or elastic (2)
if( ispec_is_acoustic( ispec_selected_source(isource) ) ) then
- idomain(isource) = 1
+ idomain(isource) = IDOMAIN_ACOUSTIC
else if( ispec_is_elastic( ispec_selected_source(isource) ) ) then
- idomain(isource) = 2
+ idomain(isource) = IDOMAIN_ELASTIC
else
idomain(isource) = 0
endif
@@ -745,31 +735,41 @@
write(IMAIN,*)
write(IMAIN,*) 'source located in slice ',islice_selected_source(isource)
write(IMAIN,*) ' in element ',ispec_selected_source(isource)
- if( idomain(isource) == 1 ) then
+
+ if( idomain(isource) == IDOMAIN_ACOUSTIC ) then
write(IMAIN,*) ' in acoustic domain'
- else if( idomain(isource) == 2 ) then
+ else if( idomain(isource) == IDOMAIN_ELASTIC ) then
write(IMAIN,*) ' in elastic domain'
else
- write(IMAIN,*) ' in unknown domain'
+ write(IMAIN,*) ' in unknown domain'
endif
write(IMAIN,*)
if(USE_FORCE_POINT_SOURCE) then
- write(IMAIN,*) ' xi coordinate of source in that element: ',nint(xi_source(isource))
- write(IMAIN,*) ' eta coordinate of source in that element: ',nint(eta_source(isource))
+ write(IMAIN,*) ' xi coordinate of source in that element: ',nint(xi_source(isource))
+ write(IMAIN,*) ' eta coordinate of source in that element: ',nint(eta_source(isource))
write(IMAIN,*) ' gamma coordinate of source in that element: ',nint(gamma_source(isource))
- write(IMAIN,*) 'nu1 = ',nu_source(1,:,isource)
- write(IMAIN,*) 'nu2 = ',nu_source(2,:,isource)
- write(IMAIN,*) 'nu3 = ',nu_source(3,:,isource)
- write(IMAIN,*) 'at (x,y,z) coordinates = ',x_found_source(isource),y_found_source(isource),z_found_source(isource)
+ write(IMAIN,*) ' nu1 = ',nu_source(1,:,isource)
+ write(IMAIN,*) ' nu2 = ',nu_source(2,:,isource)
+ write(IMAIN,*) ' nu3 = ',nu_source(3,:,isource)
+ write(IMAIN,*) ' at (x,y,z) coordinates = ',x_found_source(isource),y_found_source(isource),z_found_source(isource)
+
+ ! prints frequency content for point forces
+ f0 = hdur(isource)
+ t0_ricker = 1.2d0/f0
+ write(IMAIN,*) ' using a source of dominant frequency ',f0
+ write(IMAIN,*) ' lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+ write(IMAIN,*) ' lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+ write(IMAIN,*) ' t0 = ',t0_ricker,'t_cmt = ',t_cmt(isource)
+
else
- write(IMAIN,*) ' xi coordinate of source in that element: ',xi_source(isource)
+ write(IMAIN,*) ' xi coordinate of source in that element: ',xi_source(isource)
write(IMAIN,*) ' eta coordinate of source in that element: ',eta_source(isource)
- write(IMAIN,*) 'gamma coordinate of source in that element: ',gamma_source(isource)
+ write(IMAIN,*) ' gamma coordinate of source in that element: ',gamma_source(isource)
endif
! add message if source is a Heaviside
- if(hdur(isource) < 5.*DT) then
+ if(hdur(isource) <= 5.*DT) then
write(IMAIN,*)
write(IMAIN,*) 'Source time function is a Heaviside, convolve later'
write(IMAIN,*)
@@ -831,7 +831,7 @@
endif
! checks CMTSOLUTION format for acoustic case
- if( idomain(isource) == 1 ) then
+ if( idomain(isource) == IDOMAIN_ACOUSTIC ) then
if( Mxx(isource) /= Myy(isource) .or. Myy(isource) /= Mzz(isource) .or. &
Mxy(isource) > TINYVAL .or. Mxz(isource) > TINYVAL .or. Myz(isource) > TINYVAL ) then
write(IMAIN,*)
@@ -845,6 +845,12 @@
endif
endif
+ ! checks source domain
+ if( idomain(isource) /= IDOMAIN_ACOUSTIC .and. idomain(isource) /= IDOMAIN_ELASTIC ) then
+ ! only acoustic/elastic domain implement yet
+ call exit_MPI(myrank,'source located in unknown domain')
+ endif
+
! end of loop on all the sources
enddo
Modified: seismo/3D/SPECFEM3D/trunk/memory_eval.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/memory_eval.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/memory_eval.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -28,63 +28,120 @@
! compute the approximate amount of static memory needed to run the solver
- subroutine memory_eval(NSPEC_AB,NGLOB_AB,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh,static_memory_size)
+ subroutine memory_eval(NSPEC_AB,NGLOB_AB,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh,&
+ OCEANS,static_memory_size)
+ use create_regions_mesh_ext_par,only: NSPEC_ANISO,ispec_is_acoustic,ispec_is_elastic
+
implicit none
include "constants.h"
-! input
-! logical, intent(in) :: ATTENUATION
+ ! input
integer, intent(in) :: NSPEC_AB,NGLOB_AB
integer, intent(in) :: max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh
-
-! output
+ logical, intent(in) :: OCEANS
+ ! output
double precision, intent(out) :: static_memory_size
+ ! local parameters
+ logical :: ACOUSTIC_SIMULATION,ELASTIC_SIMULATION
-
static_memory_size = 0.d0
-! add size of each set of static arrays multiplied by the number of such arrays
+! add size of each set of arrays multiplied by the number of such arrays
-! ibool,idoubling
- static_memory_size = static_memory_size + 2.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_AB*dble(SIZE_INTEGER)
+ ! see: initialize_simulation.f90
+ ! ibool
+ static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_AB*dble(SIZE_INTEGER)
-! xix,xiy,xiz,
-! etax,etay,etaz,
-! gammax,gammay,gammaz,jacobian
-! kappavstore,muvstore
-! flag_sediments,rho_vp,rho_vs
- static_memory_size = static_memory_size + 15.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_AB*dble(CUSTOM_REAL)
+ ! xix,xiy,xiz,
+ ! etax,etay,etaz,
+ ! gammax,gammay,gammaz,jacobian
+ static_memory_size = static_memory_size + 10.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_AB*dble(CUSTOM_REAL)
-! xstore,ystore,zstore,rmass,rmass_ocean_load
- static_memory_size = static_memory_size + 5.d0*NGLOB_AB*dble(CUSTOM_REAL)
+ ! xstore,ystore,zstore
+ static_memory_size = static_memory_size + 3.d0*NGLOB_AB*dble(CUSTOM_REAL)
-! updated_dof_ocean_load,iglob_is_inner_ext_mesh
- static_memory_size = static_memory_size + 2.d0*NGLOB_AB*dble(SIZE_LOGICAL)
+ ! kappastore,mustore
+ static_memory_size = static_memory_size + 2.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_AB*dble(CUSTOM_REAL)
+
+ ! ispec_acoustic,ispec_elastic,ispec_is_poroelastic (logical)
+ static_memory_size = static_memory_size + 3.d0*NSPEC_AB*dble(SIZE_LOGICAL)
-! ispec_is_inner_ext_mesh
- static_memory_size = static_memory_size + NSPEC_AB*dble(SIZE_LOGICAL)
+ ! see: read_mesh_databases.f90
+ ! acoustic arrays
+ call any_all_l( ANY(ispec_is_acoustic), ACOUSTIC_SIMULATION )
+ if( ACOUSTIC_SIMULATION ) then
+ ! potential_acoustic, potentical_dot_acoustic, potential_dot_dot_acoustic
+ static_memory_size = static_memory_size + 3.d0*NGLOB_AB*dble(CUSTOM_REAL)
+ ! rmass_acoustic
+ static_memory_size = static_memory_size + NGLOB_AB*dble(CUSTOM_REAL)
+ ! rhostore
+ static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_AB*dble(CUSTOM_REAL)
+ endif
-! displ,veloc,accel
- static_memory_size = static_memory_size + 3.d0*dble(NDIM)*NGLOB_AB*dble(CUSTOM_REAL)
+ ! elastic arrays
+ call any_all_l( ANY(ispec_is_elastic), ELASTIC_SIMULATION )
+ if( ELASTIC_SIMULATION ) then
+ ! displacement,velocity,acceleration
+ static_memory_size = static_memory_size + 3.d0*dble(NDIM)*NGLOB_AB*dble(CUSTOM_REAL)
+
+ ! rmass
+ static_memory_size = static_memory_size + NGLOB_AB*dble(CUSTOM_REAL)
-! my_neighbours_ext_mesh,nibool_interfaces_ext_mesh
+ ! rho_vp,rho_vs
+ static_memory_size = static_memory_size + 2.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_AB*dble(CUSTOM_REAL)
+
+ ! iflag_attenaution_store
+ static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_AB*dble(SIZE_INTEGER)
+
+ ! c11store,...c66store
+ static_memory_size = static_memory_size + 21.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_ANISO*dble(CUSTOM_REAL)
+
+ if (OCEANS ) then
+ ! rmass_ocean_load
+ static_memory_size = static_memory_size + NGLOB_AB*dble(CUSTOM_REAL)
+ ! updated_dof_ocean_load
+ static_memory_size = static_memory_size + NGLOB_AB*dble(SIZE_LOGICAL)
+ endif
+ endif
+
+ ! skipping boundary surfaces
+ ! skipping free surfaces
+ ! skipping acoustic-elastic coupling surfaces
+
+ ! MPI interfaces
+ ! my_neighbours_ext_mesh,nibool_interfaces_ext_mesh
static_memory_size = static_memory_size + 2.d0*num_interfaces_ext_mesh*dble(SIZE_INTEGER)
-! ibool_interfaces_ext_mesh
- static_memory_size = static_memory_size + max_nibool_interfaces_ext_mesh*num_interfaces_ext_mesh*dble(SIZE_INTEGER)
+ ! ibool_interfaces_ext_mesh
+ static_memory_size = static_memory_size + max_nibool_interfaces_ext_mesh*num_interfaces_ext_mesh*dble(SIZE_INTEGER)
-! buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh
- static_memory_size = static_memory_size + 2.d0*dble(NDIM)*max_nibool_interfaces_ext_mesh*num_interfaces_ext_mesh*dble(CUSTOM_REAL)
+ ! MPI communications
+ ! buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh
+ static_memory_size = static_memory_size + 2.d0*dble(NDIM)*max_nibool_interfaces_ext_mesh*num_interfaces_ext_mesh*dble(CUSTOM_REAL)
-! buffer_send_scalar_ext_mesh,buffer_recv_scalar_ext_mesh
- static_memory_size = static_memory_size + 2.d0*max_nibool_interfaces_ext_mesh*num_interfaces_ext_mesh*dble(CUSTOM_REAL)
+ ! buffer_send_scalar_ext_mesh,buffer_recv_scalar_ext_mesh
+ static_memory_size = static_memory_size + 2.d0*max_nibool_interfaces_ext_mesh*num_interfaces_ext_mesh*dble(CUSTOM_REAL)
-! request_send_vector_ext_mesh,request_recv_vector_ext_mesh,request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh
- static_memory_size = static_memory_size + 4.d0*num_interfaces_ext_mesh*dble(SIZE_INTEGER)
+ ! request_send_vector_ext_mesh,request_recv_vector_ext_mesh,request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh
+ static_memory_size = static_memory_size + 4.d0*num_interfaces_ext_mesh*dble(SIZE_INTEGER)
+ ! ispec_is_inner
+ static_memory_size = static_memory_size + NSPEC_AB*dble(SIZE_LOGICAL)
+ ! skipping phase_ispec_inner_acoustic
+ ! skipping phase_ispec_inner_elastic
+
+ ! see: prepare_timerun.f90
+ ! skipping attenuation R_xx,..R_yz and epsilondev_xx,...epsilondev_yz no information yet about NSPEC_ATTENUATION_AB
+
+ ! note: no adjoint array evaluation, since it depends on SIMULATION_TYPE which can vary for each run
+ ! and is undependant of mesh databases
+
+ ! note: no dyamic arrays like for seismograms, receivers and sources,
+ ! since it depends on number of timesteps and number of stations etc., which can also vary for each run
+
end subroutine memory_eval
!
Modified: seismo/3D/SPECFEM3D/trunk/model_tomography.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/model_tomography.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/model_tomography.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -36,7 +36,8 @@
! for external tomography....
! (regular spaced, xyz-block file in ascii)
- character (len=80) :: TOMO_FILENAME = 'DATA/veryfast_tomography_abruzzo_complete.xyz'
+ !character (len=80) :: TOMO_FILENAME = 'DATA/veryfast_tomography_abruzzo_complete.xyz'
+ character (len=80) :: TOMO_FILENAME = 'DATA/tomography_model.xyz'
! model dimensions
double precision :: ORIG_X,ORIG_Y,ORIG_Z
@@ -144,6 +145,12 @@
enddo
close(27)
+
+ if( myrank == 0 ) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'tomography model: ',trim(TOMO_FILENAME)
+ write(IMAIN,*)
+ endif
end subroutine read_model_tomography
Modified: seismo/3D/SPECFEM3D/trunk/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/read_mesh_databases.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/read_mesh_databases.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -36,6 +36,7 @@
integer :: i,j,k,ispec,iglob
integer :: iinterface,ier
real(kind=CUSTOM_REAL):: minl,maxl,min_all,max_all
+ logical, dimension(:), allocatable :: iglob_is_inner
! start reading the databasesa
@@ -129,6 +130,8 @@
read(27) rmass
if( OCEANS ) then
+ ! ocean mass matrix
+ allocate(rmass_ocean_load(NGLOB_AB))
read(27) rmass_ocean_load
endif
!pll
Modified: seismo/3D/SPECFEM3D/trunk/serial.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/serial.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/serial.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -466,6 +466,8 @@
implicit none
+ include "constants.h"
+
integer sendcount, dest, sendtag, req
real(kind=CUSTOM_REAL), dimension(sendcount) :: sendbuf
@@ -481,6 +483,7 @@
implicit none
+ include "constants.h"
integer recvcount, dest, recvtag, req
real(kind=CUSTOM_REAL), dimension(recvcount) :: recvbuf
@@ -544,6 +547,8 @@
subroutine recvv_cr(recvbuf, recvcount, dest, recvtag )
implicit none
+
+ include "constants.h"
integer recvcount,dest,recvtag
real(kind=CUSTOM_REAL),dimension(recvcount) :: recvbuf
@@ -575,7 +580,9 @@
subroutine sendv_cr(sendbuf, sendcount, dest, sendtag)
implicit none
-
+
+ include "constants.h"
+
integer sendcount,dest,sendtag
real(kind=CUSTOM_REAL),dimension(sendcount) :: sendbuf
Modified: seismo/3D/SPECFEM3D/trunk/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/specfem3D_par.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/specfem3D_par.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -77,7 +77,6 @@
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: kappastore,mustore
! additional mass matrix for ocean load
-! ocean load mass matrix is always allocated statically even if no oceans
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_ocean_load
! time scheme
@@ -187,7 +186,6 @@
! MPI partition surfaces
logical, dimension(:), allocatable :: ispec_is_inner
- logical, dimension(:), allocatable :: iglob_is_inner
! maximum of the norm of the displacement
real(kind=CUSTOM_REAL) Usolidnorm,Usolidnorm_all
Modified: seismo/3D/SPECFEM3D/trunk/write_VTK_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/write_VTK_data.f90 2010-08-16 15:50:15 UTC (rev 17092)
+++ seismo/3D/SPECFEM3D/trunk/write_VTK_data.f90 2010-08-16 20:37:29 UTC (rev 17093)
@@ -49,8 +49,9 @@
character(len=256) prname_file
! write source and receiver VTK files for Paraview
- write(IMAIN,*) ' vtk file: '
- write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
+ !debug
+ !write(IMAIN,*) ' vtk file: '
+ !write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
@@ -120,8 +121,9 @@
integer :: ispec,i,j,k,ier,iglob
! write source and receiver VTK files for Paraview
- write(IMAIN,*) ' vtk file: '
- write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
+ !debug
+ !write(IMAIN,*) ' vtk file: '
+ !write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
@@ -210,8 +212,9 @@
integer :: ispec,i,j,k,ier,iglob
! write source and receiver VTK files for Paraview
- write(IMAIN,*) ' vtk file: '
- write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
+ !debug
+ !write(IMAIN,*) ' vtk file: '
+ !write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
@@ -297,8 +300,9 @@
integer :: i,iglob
! write source and receiver VTK files for Paraview
- write(IMAIN,*) ' vtk file: '
- write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
+ !debug
+ !write(IMAIN,*) ' vtk file: '
+ !write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
@@ -352,8 +356,9 @@
character(len=256) prname_file
! write source and receiver VTK files for Paraview
- write(IMAIN,*) ' vtk file: '
- write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
+ !debug
+ !write(IMAIN,*) ' vtk file: '
+ !write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
More information about the CIG-COMMITS
mailing list