[cig-commits] r18379 - in seismo/3D/FAULT_SOURCE: branches/new_fault_db branches/new_fault_db/CUBIT branches/new_fault_db/decompose_mesh_SCOTCH branches/new_fault_db/decompose_mesh_SCOTCH/devel branches/new_fault_db/run branches/new_fault_db/run/bigstar branches/new_fault_db/run/brutus branches/new_fault_db/src trunk/EXAMPLES/tpv5/DATA/FAULT
percygalvez at geodynamics.org
percygalvez at geodynamics.org
Sun May 15 19:06:15 PDT 2011
Author: percygalvez
Date: 2011-05-15 19:06:14 -0700 (Sun, 15 May 2011)
New Revision: 18379
Added:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/Makefile_bigstar
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/devel/
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/devel/decompose_mesh_SCOTCH.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/devel/fault_scotch.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/lap_to_brutus
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/bigstar/
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/bigstar/run.xgenerate_databases
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/bigstar/run.xspecfem3d
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/go_generate_databases_lsf.bash
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/go_specfem3D_lsf.bash
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/run_generate_database.bash
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/run_specfem3d.bash
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/go_generate_databases_lsf.bash
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/go_specfem3D_lsf.bash
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_generate_database.bash
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_lsf.bash
seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_specfem3d.bash
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/Makefile
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/compile_all.bash
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/tpv5/DATA/FAULT/Par_file_faults.in
Log:
new data base tpv5 test working
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py 2011-05-16 00:11:00 UTC (rev 18378)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py 2011-05-16 02:06:14 UTC (rev 18379)
@@ -19,15 +19,15 @@
ybulk=[]
zbulk=[]
-xbulk.append(-18*km) #x1
-xbulk.append(18*km) #x2
-xbulk.append(18*km) #x3
-xbulk.append(-18*km) #x4
+xbulk.append(-21*km) #x1
+xbulk.append(21*km) #x2
+xbulk.append(21*km) #x3
+xbulk.append(-21*km) #x4
-ybulk.append(-18*km) #y1
-ybulk.append(-18*km) #y2
-ybulk.append(18*km) #y3
-ybulk.append(18*km) #y4
+ybulk.append(-21*km) #y1
+ybulk.append(-21*km) #y2
+ybulk.append(21*km) #y3
+ybulk.append(21*km) #y4
zbulk=[z_surf]*4
@@ -85,7 +85,7 @@
cubit.cmd("merge all")
cubit.cmd("surface 1 size "+str(elementsize))
cubit.cmd("volume 1 size "+str(elementsize))
-
+cubit.cmd("surface 1 scheme pave")
cubit.cmd("mesh surface 1")
cubit.cmd("mesh volume 1")
#cubit.cmd("unmerge surface 2 3")
@@ -128,7 +128,7 @@
for f in faces:
if dic_quads_fault_u.has_key(f):
nodes=cubit.get_connectivity('Face',f)
- print 'h,fault nodes side up :',h,nodes[0],nodes[1],nodes[2],nodes[3]
+# print 'h,fault nodes side up :',h,nodes[0],nodes[1],nodes[2],nodes[3]
txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
nodes[1],nodes[2],nodes[3])
fault_file.write(txt)
@@ -139,7 +139,7 @@
for f in faces:
if dic_quads_fault_d.has_key(f):
nodes=cubit.get_connectivity('Face',f)
- print 'h,fault nodes side down :',h,nodes[0],nodes[1],nodes[2],nodes[3]
+# print 'h,fault nodes side down :',h,nodes[0],nodes[1],nodes[2],nodes[3]
txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
nodes[1],nodes[2],nodes[3])
fault_file.write(txt)
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/Makefile
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/Makefile 2011-05-16 00:11:00 UTC (rev 18378)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/Makefile 2011-05-16 02:06:14 UTC (rev 18379)
@@ -6,7 +6,7 @@
F90 = gfortran # use -Wall
## modify to match your library paths
-SCOTCH_LIBS = -L/home/galvez/scotch/scotch_5.1/lib/ -lscotch -lscotcherr
+SCOTCH_LIBS = -L/cluster/home/erdw/gpercy/scotch/scotch_5.1/lib/ -lscotch -lscotcherr
#(SCOTCH library address .Change it in case
#your SCOTH is install in another address)
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/Makefile_bigstar
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/Makefile_bigstar (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/Makefile_bigstar 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,39 @@
+# Makefile
+
+#############################################################
+## Compiler :
+#F90 = ifort # use -warn
+F90 = gfortran # use -Wall
+
+## modify to match your library paths
+SCOTCH_LIBS = -L/home/galvez/scotch/scotch_5.1/lib/ -lscotch -lscotcherr
+#(SCOTCH library address .Change it in case
+#your SCOTH is install in another address)
+
+#############################################################
+
+LIBS = fault_scotch.o\
+ part_decompose_mesh_SCOTCH.o \
+ decompose_mesh_SCOTCH.o \
+ program_decompose_mesh_SCOTCH.o
+
+# targets
+all: xdecompose_mesh_SCOTCH
+
+xdecompose_mesh_SCOTCH: $(LIBS)
+ ${F90} -o xdecompose_mesh_SCOTCH $(LIBS) $(SCOTCH_LIBS)
+
+fault_scotch.o : fault_scotch.f90
+ ${F90} -c fault_scotch.f90
+
+part_decompose_mesh_SCOTCH.o: part_decompose_mesh_SCOTCH.f90
+ ${F90} -c part_decompose_mesh_SCOTCH.f90
+
+decompose_mesh_SCOTCH.o: decompose_mesh_SCOTCH.f90 part_decompose_mesh_SCOTCH.f90
+ ${F90} -c decompose_mesh_SCOTCH.f90
+
+program_decompose_mesh_SCOTCH.o: program_decompose_mesh_SCOTCH.f90
+ ${F90} -c program_decompose_mesh_SCOTCH.f90
+
+clean:
+ rm -f *.o *.mod a.out xdecompose_mesh_SCOTCH
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/compile_all.bash
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/compile_all.bash 2011-05-16 00:11:00 UTC (rev 18378)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/compile_all.bash 2011-05-16 02:06:14 UTC (rev 18379)
@@ -2,7 +2,5 @@
make clean
make
# ./xdecompose_mesh_SCOTCH 4 MESH/OUTPUT_FILES ../DATABASES_MPI/
-#npart=`grep nparts constants_decompose_mesh_SCOTCH.h | cut -d = -f 2`
-npart = 4
echo "./xdecompose_mesh_SCOTCH $npart $input_dir $ouput_dir "
-./xdecompose_mesh_SCOTCH $npart ../CUBIT/MESH/ ../DATABASES_MPI/
+./xdecompose_mesh_SCOTCH 100 ~/FAULT_SOURCE/CUBIT/MESH /cluster/work/erdw/gpercy
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/devel/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/devel/decompose_mesh_SCOTCH.f90 (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/devel/decompose_mesh_SCOTCH.f90 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,699 @@
+!program pre_meshfem3D
+
+module decompose_mesh_SCOTCH
+
+ use part_decompose_mesh_SCOTCH
+ use fault_scotch
+
+ implicit none
+
+ include './scotchf.h'
+
+! number of partitions
+ integer :: nparts ! e.g. 4 for partitioning for 4 CPUs or 4 processes
+
+! mesh arrays
+ integer(long) :: nspec
+ integer, dimension(:,:), allocatable :: elmnts
+ integer, dimension(:,:), allocatable :: mat
+ integer, dimension(:), allocatable :: part
+
+ integer :: nnodes
+ double precision, dimension(:,:), allocatable :: nodes_coords
+
+ integer, dimension(:), allocatable :: xadj
+ integer, dimension(:), allocatable :: adjncy
+ integer, dimension(:), allocatable :: nnodes_elmnts
+ integer, dimension(:), allocatable :: nodes_elmnts
+ integer, dimension(:), allocatable :: elmnts_load
+
+ integer, dimension(:), pointer :: glob2loc_elmnts
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+
+ integer, dimension(:), pointer :: tab_size_interfaces, tab_interfaces
+ integer, dimension(:), allocatable :: my_interfaces
+ integer, dimension(:), allocatable :: my_nb_interfaces
+ integer :: ninterfaces
+ integer :: my_ninterface
+
+ integer(long) :: nsize ! Max number of elements that contain the same node.
+ integer :: nb_edges
+
+ integer :: ispec, inode
+ integer :: ngnod
+ integer :: max_neighbour ! Real maximum number of neighbours per element
+ integer(long) :: sup_neighbour ! Majoration of the maximum number of neighbours per element
+
+ integer :: ipart, nnodes_loc, nspec_loc
+ integer :: num_elmnt, num_node, num_mat
+
+ ! boundaries
+ integer :: ispec2D
+ integer :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, nspec2D_bottom, nspec2D_top
+ integer, dimension(:), allocatable :: ibelm_xmin, ibelm_xmax, ibelm_ymin, ibelm_ymax, ibelm_bottom, ibelm_top
+ integer, dimension(:,:), allocatable :: nodes_ibelm_xmin, nodes_ibelm_xmax, nodes_ibelm_ymin
+ integer, dimension(:,:), allocatable :: nodes_ibelm_ymax, nodes_ibelm_bottom, nodes_ibelm_top
+
+ ! moho surface (optional)
+ integer :: nspec2D_moho
+ integer, dimension(:), allocatable :: ibelm_moho
+ integer, dimension(:,:), allocatable :: nodes_ibelm_moho
+
+ character(len=256) :: prname
+
+ logical, dimension(:), allocatable :: mask_nodes_elmnts
+ integer, dimension(:), allocatable :: used_nodes_elmnts
+
+ double precision, dimension(SCOTCH_GRAPHDIM) :: scotchgraph
+ double precision, dimension(SCOTCH_STRATDIM) :: scotchstrat
+ character(len=256), parameter :: scotch_strategy='b{job=t,map=t,poli=S,sep=h{pass=30}}'
+ integer :: ierr,idummy
+
+ !pll
+ double precision , dimension(:,:), allocatable :: mat_prop
+ integer :: count_def_mat,count_undef_mat,imat
+ character (len=30), dimension(:,:), allocatable :: undef_mat_prop
+
+! default mesh file directory
+ character(len=256) :: localpath_name
+ character(len=256) :: outputpath_name
+
+ integer :: q_flag,aniso_flag,idomain_id
+ double precision :: vp,vs,rho
+
+ contains
+
+ !----------------------------------------------------------------------------------------------
+ ! reads in mesh files
+ !----------------------------------------------------------------------------------------------
+ subroutine read_mesh_files
+
+ ! sets number of nodes per element
+ ngnod = esize
+
+ ! reads node coordinates
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/nodes_coords_file',&
+ status='old', form='formatted', iostat = ierr)
+ if( ierr /= 0 ) then
+ print*,'could not open file:',localpath_name(1:len_trim(localpath_name))//'/nodes_coords_file'
+ stop 'error file open'
+ endif
+ read(98,*) nnodes
+ allocate(nodes_coords(3,nnodes))
+ do inode = 1, nnodes
+ ! format: #id_node #x_coordinate #y_coordinate #z_coordinate
+ read(98,*) num_node, nodes_coords(1,num_node), nodes_coords(2,num_node), nodes_coords(3,num_node)
+ !if(num_node /= inode) stop "ERROR : Invalid nodes_coords file."
+ end do
+ close(98)
+ print*, 'total number of nodes: '
+ print*, ' nnodes = ', nnodes
+
+ ! reads mesh elements indexing
+ !(CUBIT calls this the connectivity, guess in the sense that it connects with the points index in
+ ! the global coordinate file "nodes_coords_file"; it doesn't tell you which point is connected with others)
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/mesh_file', &
+ status='old', form='formatted')
+ read(98,*) nspec
+ allocate(elmnts(esize,nspec))
+ do ispec = 1, nspec
+ ! format: # element_id #id_node1 ... #id_node8
+
+ ! note: be aware that here we can have different node ordering for a cube element;
+ ! the ordering from Cubit files might not be consistent for multiple volumes, or uneven, unstructured grids
+ !
+ ! guess here it assumes that spectral elements ordering is like first at the bottom of the element, anticlock-wise, i.e.
+ ! point 1 = (0,0,0), point 2 = (0,1,0), point 3 = (1,1,0), point 4 = (1,0,0)
+ ! then top (positive z-direction) of element
+ ! point 5 = (0,0,1), point 6 = (0,1,1), point 7 = (1,1,1), point 8 = (1,0,1)
+
+ !read(98,*) num_elmnt, elmnts(5,num_elmnt), elmnts(1,num_elmnt),elmnts(4,num_elmnt), elmnts(8,num_elmnt), &
+ ! elmnts(6,num_elmnt), elmnts(2,num_elmnt), elmnts(3,num_elmnt), elmnts(7,num_elmnt)
+
+ read(98,*) num_elmnt, elmnts(1,num_elmnt), elmnts(2,num_elmnt),elmnts(3,num_elmnt), elmnts(4,num_elmnt), &
+ elmnts(5,num_elmnt), elmnts(6,num_elmnt), elmnts(7,num_elmnt), elmnts(8,num_elmnt)
+
+ if((num_elmnt > nspec) .or. (num_elmnt < 1) ) stop "ERROR : Invalid mesh file."
+
+
+ !outputs info for each element to see ordering
+ !print*,'ispec: ',ispec
+ !print*,' ',num_elmnt, elmnts(5,num_elmnt), elmnts(1,num_elmnt),elmnts(4,num_elmnt), elmnts(8,num_elmnt), &
+ ! elmnts(6,num_elmnt), elmnts(2,num_elmnt), elmnts(3,num_elmnt), elmnts(7,num_elmnt)
+ !print*,'elem:',num_elmnt
+ !do i=1,8
+ ! print*,' i ',i,'val :',elmnts(i,num_elmnt),&
+ ! nodes_coords(1,elmnts(i,num_elmnt)),nodes_coords(2,elmnts(i,num_elmnt)),nodes_coords(3,elmnts(i,num_elmnt))
+ !enddo
+ !print*
+
+ end do
+ close(98)
+ print*, 'total number of spectral elements:'
+ print*, ' nspec = ', nspec
+
+ ! reads material associations
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/materials_file', &
+ status='old', form='formatted')
+ allocate(mat(2,nspec))
+ do ispec = 1, nspec
+ ! format: # id_element #flag
+ ! note: be aware that elements may not be sorted in materials_file
+ read(98,*) num_mat,mat(1,num_mat) !mat(1,ispec)!, mat(2,ispec)
+ if((num_mat > nspec) .or. (num_mat < 1) ) stop "ERROR : Invalid mat file."
+ end do
+ close(98)
+
+ ! TODO:
+ ! must be changed, if mat(1,i) < 0 1 == interface , 2 == tomography
+ mat(2,:) = 1
+
+ ! reads material definitions
+ !
+ ! note: format of nummaterial_velocity_file must be
+ !
+ ! #(1)material_domain_id #(2)material_id #(3)rho #(4)vp #(5)vs #(6)Q_flag #(7)anisotropy_flag
+ !
+ ! where
+ ! material_domain_id : 1=acoustic / 2=elastic / 3=poroelastic
+ ! material_id : number of material/volume
+ ! rho : density
+ ! vp : P-velocity
+ ! vs : S-velocity
+ ! Q_flag : 0=no attenuation/1=IATTENUATION_SEDIMENTS_40, 2=..., 13=IATTENUATION_BEDROCK
+ ! anisotropy_flag : 0=no anisotropy/ 1,2,.. check with implementation in aniso_model.f90
+ count_def_mat = 0
+ count_undef_mat = 0
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/nummaterial_velocity_file',&
+ status='old', form='formatted')
+ ! note: format #material_domain_id #material_id #...
+ read(98,*,iostat=ierr) idummy,num_mat
+ print *,'materials:'
+ ! counts materials (defined/undefined)
+ do while (ierr == 0)
+ print*, ' num_mat = ',num_mat
+ if(num_mat /= -1) then
+ count_def_mat = count_def_mat + 1
+ else
+ count_undef_mat = count_undef_mat + 1
+ end if
+ read(98,*,iostat=ierr) idummy,num_mat
+ end do
+ close(98)
+ print*, ' defined = ',count_def_mat, 'undefined = ',count_undef_mat
+ ! check with material flags
+ if( count_def_mat > 0 .and. maxval(mat(1,:)) > count_def_mat ) then
+ print*,'error material definitions:'
+ print*,' materials associated in materials_file:',maxval(mat(1,:))
+ print*,' bigger than defined materials in nummaterial_velocity_file:',count_def_mat
+ stop 'error materials'
+ endif
+ allocate(mat_prop(6,count_def_mat))
+ allocate(undef_mat_prop(6,count_undef_mat))
+ ! reads in defined material properties
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/nummaterial_velocity_file', &
+ status='old', form='formatted')
+ do imat=1,count_def_mat
+ ! material definitions
+ !
+ ! format: note that we save the arguments in a slightly different order in mat_prop(:,:)
+ ! #(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
+ !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
+ mat_prop(2,num_mat) = vp
+ 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."
+
+ !checks attenuation flag with integer range as defined in constants.h like IATTENUATION_SEDIMENTS_40, ....
+ if( int(mat_prop(4,num_mat)) > 13 ) then
+ stop 'wrong attenuation flag in mesh: too large, not supported yet - check with constants.h'
+ endif
+ end do
+ ! 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),&
+ undef_mat_prop(3,imat),undef_mat_prop(4,imat),undef_mat_prop(5,imat)
+ end do
+ close(98)
+
+ ! 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)
+ if( ierr /= 0 ) then
+ nspec2D_xmin = 0
+ else
+ read(98,*) nspec2D_xmin
+ endif
+ allocate(ibelm_xmin(nspec2D_xmin))
+ allocate(nodes_ibelm_xmin(4,nspec2D_xmin))
+ do ispec2D = 1,nspec2D_xmin
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ ! note: ordering for CUBIT seems such that the normal of the face points outward of the element the face belongs to;
+ ! in other words, nodes are in increasing order such that when looking from within the element outwards,
+ ! they are ordered clockwise
+ !
+ ! doesn't necessarily have to start on top-rear, then bottom-rear, bottom-front, and finally top-front i.e.:
+ ! point 1 = (0,1,1), point 2 = (0,1,0), point 3 = (0,0,0), point 4 = (0,0,1)
+ read(98,*) ibelm_xmin(ispec2D), nodes_ibelm_xmin(1,ispec2D), nodes_ibelm_xmin(2,ispec2D), &
+ nodes_ibelm_xmin(3,ispec2D), nodes_ibelm_xmin(4,ispec2D)
+
+ !outputs info for each element for check of ordering
+ !print*,'ispec2d:',ispec2d
+ !print*,' xmin:', ibelm_xmin(ispec2D), nodes_ibelm_xmin(1,ispec2D), nodes_ibelm_xmin(2,ispec2D), &
+ ! nodes_ibelm_xmin(3,ispec2D), nodes_ibelm_xmin(4,ispec2D)
+ !do i=1,4
+ ! print*,'i',i,'val:',ibelm_xmin(ispec2d),nodes_coords(1,nodes_ibelm_xmin(i,ispec2D)), &
+ ! nodes_coords(2,nodes_ibelm_xmin(i,ispec2D)),nodes_coords(3,nodes_ibelm_xmin(i,ispec2D))
+ !enddo
+ !print*
+ end do
+ close(98)
+ print*, 'absorbing boundaries:'
+ print*, ' nspec2D_xmin = ', nspec2D_xmin
+
+ ! reads in absorbing boundary files
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_xmax', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_xmax = 0
+ else
+ read(98,*) nspec2D_xmax
+ endif
+ allocate(ibelm_xmax(nspec2D_xmax))
+ allocate(nodes_ibelm_xmax(4,nspec2D_xmax))
+ do ispec2D = 1,nspec2D_xmax
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ read(98,*) ibelm_xmax(ispec2D), nodes_ibelm_xmax(1,ispec2D), nodes_ibelm_xmax(2,ispec2D), &
+ nodes_ibelm_xmax(3,ispec2D), nodes_ibelm_xmax(4,ispec2D)
+ end do
+ close(98)
+ print*, ' nspec2D_xmax = ', nspec2D_xmax
+
+ ! reads in absorbing boundary files
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_ymin', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_ymin = 0
+ else
+ read(98,*) nspec2D_ymin
+ endif
+ allocate(ibelm_ymin(nspec2D_ymin))
+ allocate(nodes_ibelm_ymin(4,nspec2D_ymin))
+ do ispec2D = 1,nspec2D_ymin
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ read(98,*) ibelm_ymin(ispec2D), nodes_ibelm_ymin(1,ispec2D), nodes_ibelm_ymin(2,ispec2D), &
+ nodes_ibelm_ymin(3,ispec2D), nodes_ibelm_ymin(4,ispec2D)
+ end do
+ close(98)
+ print*, ' nspec2D_ymin = ', nspec2D_ymin
+
+ ! reads in absorbing boundary files
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_ymax', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_ymax = 0
+ else
+ read(98,*) nspec2D_ymax
+ endif
+ allocate(ibelm_ymax(nspec2D_ymax))
+ allocate(nodes_ibelm_ymax(4,nspec2D_ymax))
+ do ispec2D = 1,nspec2D_ymax
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ read(98,*) ibelm_ymax(ispec2D), nodes_ibelm_ymax(1,ispec2D), nodes_ibelm_ymax(2,ispec2D), &
+ nodes_ibelm_ymax(3,ispec2D), nodes_ibelm_ymax(4,ispec2D)
+ end do
+ close(98)
+ print*, ' nspec2D_ymax = ', nspec2D_ymax
+
+ ! reads in absorbing boundary files
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_bottom', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_bottom = 0
+ else
+ read(98,*) nspec2D_bottom
+ endif
+ allocate(ibelm_bottom(nspec2D_bottom))
+ allocate(nodes_ibelm_bottom(4,nspec2D_bottom))
+ do ispec2D = 1,nspec2D_bottom
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ read(98,*) ibelm_bottom(ispec2D), nodes_ibelm_bottom(1,ispec2D), nodes_ibelm_bottom(2,ispec2D), &
+ nodes_ibelm_bottom(3,ispec2D), nodes_ibelm_bottom(4,ispec2D)
+ end do
+ close(98)
+ print*, ' nspec2D_bottom = ', nspec2D_bottom
+
+ ! reads in free_surface boundary files
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/free_surface_file', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_top = 0
+ else
+ read(98,*) nspec2D_top
+ endif
+ allocate(ibelm_top(nspec2D_top))
+ allocate(nodes_ibelm_top(4,nspec2D_top))
+ do ispec2D = 1,nspec2D_top
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ read(98,*) ibelm_top(ispec2D), nodes_ibelm_top(1,ispec2D), nodes_ibelm_top(2,ispec2D), &
+ nodes_ibelm_top(3,ispec2D), nodes_ibelm_top(4,ispec2D)
+ end do
+ close(98)
+ print*, ' nspec2D_top = ', nspec2D_top
+
+ ! reads in moho_surface boundary files (optional)
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/moho_surface_file', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_moho = 0
+ else
+ read(98,*) nspec2D_moho
+ endif
+ allocate(ibelm_moho(nspec2D_moho))
+ allocate(nodes_ibelm_moho(4,nspec2D_moho))
+ do ispec2D = 1,nspec2D_moho
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ read(98,*) ibelm_moho(ispec2D), nodes_ibelm_moho(1,ispec2D), nodes_ibelm_moho(2,ispec2D), &
+ nodes_ibelm_moho(3,ispec2D), nodes_ibelm_moho(4,ispec2D)
+ end do
+ close(98)
+ if( nspec2D_moho > 0 ) print*, ' nspec2D_moho = ', nspec2D_moho
+
+ call read_fault_files(localpath_name)
+ if (allocated(faults)) then
+ call save_nodes_coords(nodes_coords,nnodes)
+ call close_faults(nodes_coords,nnodes)
+ end if
+
+ end subroutine read_mesh_files
+
+ !----------------------------------------------------------------------------------------------
+ ! checks valence of nodes
+ !----------------------------------------------------------------------------------------------
+
+ subroutine check_valence
+
+ allocate(mask_nodes_elmnts(nnodes))
+ allocate(used_nodes_elmnts(nnodes))
+ mask_nodes_elmnts(:) = .false.
+ used_nodes_elmnts(:) = 0
+ do ispec = 1, nspec
+ do inode = 1, ESIZE
+ mask_nodes_elmnts(elmnts(inode,ispec)) = .true.
+ used_nodes_elmnts(elmnts(inode,ispec)) = used_nodes_elmnts(elmnts(inode,ispec)) + 1
+ enddo
+ enddo
+ print *, 'nodes valence: '
+ print *, ' min = ',minval(used_nodes_elmnts(:)),'max = ', maxval(used_nodes_elmnts(:))
+ do inode = 1, nnodes
+ if (.not. mask_nodes_elmnts(inode)) then
+ stop 'ERROR : nodes not used.'
+ endif
+ enddo
+ nsize = maxval(used_nodes_elmnts(:))
+ sup_neighbour = ngnod * nsize - (ngnod + (ngnod/2 - 1)*nfaces)
+ print*, ' nsize = ',nsize, 'sup_neighbour = ', sup_neighbour
+
+ end subroutine check_valence
+
+ !----------------------------------------------------------------------------------------------
+ ! divides model into partitions using scotch library functions
+ !----------------------------------------------------------------------------------------------
+
+ subroutine scotch_partitioning
+
+ implicit none
+
+ elmnts(:,:) = elmnts(:,:) - 1
+
+ ! determines maximum neighbors based on 1 common node
+ allocate(xadj(1:nspec+1))
+ allocate(adjncy(1:sup_neighbour*nspec))
+ allocate(nnodes_elmnts(1:nnodes))
+ allocate(nodes_elmnts(1:nsize*nnodes))
+ call mesh2dual_ncommonnodes(nspec, nnodes, nsize, sup_neighbour, elmnts, xadj, adjncy, nnodes_elmnts, &
+ nodes_elmnts, max_neighbour, 1)
+ print*, 'mesh2dual: '
+ print*, ' max_neighbour = ',max_neighbour
+
+ nb_edges = xadj(nspec+1)
+
+ ! allocates & initializes partioning of elements
+ allocate(part(1:nspec))
+ part(:) = -1
+
+ ! initializes
+ ! elements load array
+ allocate(elmnts_load(1:nspec))
+
+ ! uniform load
+ 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)
+
+ ! SCOTCH partitioning
+ call scotchfstratinit (scotchstrat(1), ierr)
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Cannot initialize 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'
+ endif
+
+ ! fills graph structure : see user manual (scotch_user5.1.pdf, page 72/73)
+ ! arguments: #(1) graph_structure #(2) baseval(either 0/1) #(3) number_of_vertices
+ ! #(4) adjacency_index_array #(5) adjacency_end_index_array (optional)
+ ! #(6) vertex_load_array (optional) #(7) vertex_label_array
+ ! #(7) number_of_arcs #(8) adjacency_array
+ ! #(9) arc_load_array (optional) #(10) ierror
+ call scotchfgraphbuild (scotchgraph (1), 0, nspec, &
+ xadj (1), xadj (1), &
+ elmnts_load (1), xadj (1), &
+ nb_edges, adjncy (1), &
+ adjncy (1), ierr)
+
+ ! w/out element load, but adjacency array
+ !call scotchfgraphbuild (scotchgraph (1), 0, nspec, &
+ ! xadj (1), xadj (1), &
+ ! xadj (1), xadj (1), &
+ ! nb_edges, adjncy (1), &
+ ! adjncy (1), ierr)
+
+
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Cannot build graph'
+ endif
+
+ call scotchfgraphcheck (scotchgraph (1), ierr)
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Invalid check'
+ endif
+
+ call scotchfgraphpart (scotchgraph (1), nparts, scotchstrat(1),part(1),ierr)
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Cannot part graph'
+ endif
+
+ call scotchfgraphexit (scotchgraph (1), ierr)
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Cannot destroy graph'
+ endif
+
+ call scotchfstratexit (scotchstrat(1), ierr)
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Cannot destroy strat'
+ endif
+
+
+ ! re-partitioning puts acoustic-elastic coupled elements into same partition
+ ! integer :: nfaces_coupled
+ ! integer, dimension(:,:), pointer :: faces_coupled
+ ! call acoustic_elastic_repartitioning (nspec, nnodes, elmnts, &
+ ! count_def_mat, mat(1,:) , mat_prop, &
+ ! sup_neighbour, nsize, &
+ ! nparts, part, nfaces_coupled, faces_coupled)
+
+ ! move all fault elements to the same partition (proc=0)
+ call fault_collect_elements(nspec,nnodes,elmnts, &
+ sup_neighbour,esize,nsize,nparts,part)
+
+ ! re-partitioning puts moho-surface coupled elements into same partition
+ call moho_surface_repartitioning (nspec, nnodes, elmnts, &
+ sup_neighbour, nsize, nparts, part, &
+ nspec2D_moho,ibelm_moho,nodes_ibelm_moho )
+
+
+ ! local number of each element for each partition
+ call Construct_glob2loc_elmnts(nspec, part, glob2loc_elmnts,nparts)
+
+
+ ! local number of each node for each partition
+ call Construct_glob2loc_nodes(nspec, nnodes,nsize, nnodes_elmnts, nodes_elmnts, part, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, nparts)
+
+
+ ! mpi interfaces
+ ! acoustic/elastic boundaries WILL BE SEPARATED into different MPI partitions
+ call Construct_interfaces(nspec, sup_neighbour, part, elmnts, &
+ xadj, adjncy, tab_interfaces, &
+ tab_size_interfaces, ninterfaces, &
+ nparts)
+
+ !or: uncomment if you want acoustic/elastic boundaries NOT to be separated into different MPI partitions
+ !call Construct_interfaces_no_ac_el_sep(nspec, sup_neighbour, part, elmnts, &
+ ! xadj, adjncy, tab_interfaces, &
+ ! tab_size_interfaces, ninterfaces, &
+ ! count_def_mat, mat_prop(3,:), mat(1,:), nparts)
+
+
+ end subroutine scotch_partitioning
+
+
+ !----------------------------------------------------------------------------------------------
+ ! writes out new Databases files for each partition
+ !----------------------------------------------------------------------------------------------
+
+ subroutine write_mesh_databases
+
+ allocate(my_interfaces(0:ninterfaces-1))
+ allocate(my_nb_interfaces(0:ninterfaces-1))
+
+ ! writes out Database file for each partition
+
+ do ipart = 0, nparts-1
+
+ ! opens output file
+ write(prname, "(i6.6,'_Database')") ipart
+ open(unit=15,file=outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname,&
+ status='unknown', action='write', form='formatted', iostat = ierr)
+ if( ierr /= 0 ) then
+ print*,'error file open:',outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname
+ print*
+ print*,'check if path exists:',outputpath_name(1:len_trim(outputpath_name))
+ stop 'error file open Database'
+ endif
+
+ ! gets number of nodes
+ call write_glob2loc_nodes_database(15, ipart, nnodes_loc, nodes_coords, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, nnodes, 1)
+
+ ! gets number of spectral elements
+ call write_partition_database(15, ipart, nspec_loc, nspec, elmnts, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part, mat, ngnod, 1)
+
+ ! writes out node coordinate locations
+ write(15,*) nnodes_loc
+
+ call write_glob2loc_nodes_database(15, ipart, nnodes_loc, nodes_coords,&
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, nnodes, 2)
+
+ call write_material_properties_database(15,count_def_mat,count_undef_mat, &
+ mat_prop, undef_mat_prop)
+
+ ! writes out spectral element indices
+ write(15,*) nspec_loc
+
+ call write_partition_database(15, ipart, nspec_loc, nspec, elmnts, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part, mat, ngnod, 2)
+
+ ! writes out absorbing/free-surface boundaries
+ call write_boundaries_database(15, ipart, nspec, nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, &
+ nspec2D_ymax, nspec2D_bottom, nspec2D_top, &
+ ibelm_xmin, ibelm_xmax, ibelm_ymin, &
+ ibelm_ymax, ibelm_bottom, ibelm_top, &
+ nodes_ibelm_xmin, nodes_ibelm_xmax, nodes_ibelm_ymin, &
+ nodes_ibelm_ymax, nodes_ibelm_bottom, nodes_ibelm_top, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part)
+
+ ! gets number of MPI interfaces
+ call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, ipart, ninterfaces, &
+ my_ninterface, my_interfaces, my_nb_interfaces, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, 1, nparts)
+
+ ! writes out MPI interfaces elements
+ write(15,*) my_ninterface, maxval(my_nb_interfaces)
+
+ call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, ipart, ninterfaces, &
+ my_ninterface, my_interfaces, my_nb_interfaces, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, 2, nparts)
+
+ ! writes out moho surface (optional)
+ call write_moho_surface_database(15, ipart, nspec, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part, &
+ nspec2D_moho,ibelm_moho,nodes_ibelm_moho)
+
+ close(15)
+
+ if (.not. allocated(faults)) cycle
+ ! write fault database
+ write(prname, "(i6.6,'_Database_fault')") ipart
+ open(unit=16,file=outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname,&
+ status='unknown', action='write', form='formatted', iostat = ierr)
+ if( ierr /= 0 ) then
+ print*,'error file open:',outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname
+ print*
+ print*,'check if path exists:',outputpath_name(1:len_trim(outputpath_name))
+ stop
+ endif
+ call write_fault_database(16, ipart, nspec, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, part)
+! close(16)
+
+ ! write nodes coordinates with fault open crack
+! write(prname, "(i6.6,'_Database_fault_nodes')") ipart
+! open(unit=16,file=outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname,&
+! status='unknown', action='write', form='formatted', iostat = ierr)
+! if( ierr /= 0 ) then
+! print*,'error file open:',outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname
+! print*
+! print*,'check if path exists:',outputpath_name(1:len_trim(outputpath_name))
+! stop
+! endif
+
+ write(16,*) nnodes_loc
+ call write_glob2loc_nodes_database(16, ipart, nnodes_loc, nodes_coords_open,&
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, nnodes, 2)
+
+ close(16)
+
+ end do
+ print*, 'partitions: '
+ print*, ' num = ',nparts
+ print*
+ print*, 'Databases files in directory: ',outputpath_name(1:len_trim(outputpath_name))
+ print*, 'finished successfully'
+ print*
+
+ end subroutine write_mesh_databases
+
+!end program pre_meshfem3D
+
+end module decompose_mesh_SCOTCH
+
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/devel/fault_scotch.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/devel/fault_scotch.f90 (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/devel/fault_scotch.f90 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,581 @@
+module fault_scotch
+
+ implicit none
+
+ private
+ type fault_type
+ private
+ integer :: nspec
+ integer, dimension(:), pointer :: ispec1, ispec2 !, iface1, iface2
+ integer, dimension(:,:), pointer :: inodes1, inodes2
+ end type fault_type
+
+ type(fault_type), allocatable, save :: faults(:)
+ double precision, dimension(:,:), allocatable, save :: nodes_coords_open
+
+
+ integer, parameter :: long = SELECTED_INT_KIND(18)
+
+ double precision, parameter :: FAULT_GAP_TOLERANCE = 1.0d0 ! JPA: are you sure 1 meter is small enough?
+ ! PGB: for a simple test is fine .For SCEC lower values.
+
+ public :: read_fault_files, fault_collect_elements, close_faults, write_fault_database, &
+ save_nodes_coords, nodes_coords_open, faults
+
+CONTAINS
+!==========================================================================================
+
+ Subroutine read_fault_files(localpath_name)
+
+ character(len=256),intent(in) :: localpath_name
+ integer :: nbfaults, iflt, ier
+
+ open(101,file='../DATA/FAULT/Par_file_faults.in',status='old',action='read',iostat=ier)
+ if (ier==0) then
+ read(101,*) nbfaults
+ else
+ nbfaults = 0
+ print *, 'Par_file.in not found: assume no faults'
+ endif
+ close(101)
+
+ if (nbfaults>0) then
+ allocate(faults(nbfaults))
+ do iflt = 1 , nbfaults
+ call read_single_fault_file(faults(iflt),iflt,localpath_name)
+ enddo
+ endif
+
+ end subroutine read_fault_files
+
+
+!---------------------------------------------------------------------------------------------------
+
+ Subroutine read_single_fault_file(f,ifault,localpath_name)
+
+ type(fault_type), intent(inout) :: f
+ character(len=256),intent(in) :: localpath_name
+
+ character(len=256) :: filename
+ integer,intent(in) :: ifault
+ character(len=5) :: NTchar
+ integer :: e,ier,nspec_side1,nspec_side2
+
+ write(NTchar,'(I5)') ifault
+ NTchar = adjustl(NTchar)
+
+ filename = localpath_name(1:len_trim(localpath_name))//'/fault_file_'//&
+ NTchar(1:len_trim(NTchar))//'.dat'
+ filename = adjustl(filename)
+ ! reads fault elements and nodes
+ open(unit=101, file=filename, status='old', form='formatted', iostat = ier)
+ if( ier /= 0 ) then
+ write(6,*) 'Fatal error: file '//filename//' not found'
+ write(6,*) 'Abort'
+ stop
+ endif
+
+ read(101,*) nspec_side1,nspec_side2
+ if (nspec_side1 /= nspec_side2) stop 'Number of fault nodes at do not match.'
+ f%nspec = nspec_side1
+ allocate(f%ispec1(f%nspec))
+ allocate(f%ispec2(f%nspec))
+ allocate(f%inodes1(4,f%nspec))
+ allocate(f%inodes2(4,f%nspec))
+ ! format:
+ ! number of elements side at side 1 and side 2
+ ! #id_(element containing the face) #id_node1_face .. #id_node4_face
+ ! First for all faces on side 1, then side 2
+ ! Note: element ids start at 1, not 0. Check that in cubit2specfem3d.py
+ do e=1,f%nspec
+ read(101,*) f%ispec1(e),f%inodes1(:,e)
+! print*, f%ispec1(e),f%inodes1(:,e) !TEST
+ enddo
+ do e=1,f%nspec
+ read(101,*) f%ispec2(e),f%inodes2(:,e)
+! print*, f%ispec2(e),f%inodes2(:,e) !TEST
+ enddo
+ ! If we ever figure out how to export "ifaces" from CUBIT:
+ !allocate(f%iface1(f%nspec))
+ !allocate(f%iface2(f%nspec))
+ !do e=1,f%nspec
+ ! read(101,*) f%ispec1(e),f%ispec2(e),f%iface1(e),f%iface2(e)
+ !enddo
+
+ close(101)
+
+ end Subroutine read_single_fault_file
+
+
+! ---------------------------------------------------------------------------------------------------
+! Saving nodes_coords to be used in SESAME for ibool_fault_side1 and side2
+ subroutine save_nodes_coords(nodes_coords,nnodes)
+
+ integer, intent(in) :: nnodes
+ double precision, dimension(3,nnodes), intent(in) :: nodes_coords
+
+ allocate(nodes_coords_open(3,nnodes))
+ nodes_coords_open = nodes_coords
+
+ end subroutine save_nodes_coords
+
+! ---------------------------------------------------------------------------------------------------
+
+ subroutine close_faults(nodes_coords,nnodes)
+
+ integer ,intent(in) :: nnodes
+ double precision, dimension(3,nnodes), intent(inout) :: nodes_coords
+
+ integer :: iflt
+
+ do iflt=1,size(faults)
+ call close_fault_single(faults(iflt),nodes_coords,nnodes)
+ enddo
+
+ end subroutine close_faults
+
+! ---------------------------------------------------------------------------------------------------
+!jpa: to do this much faster:
+! 1. create a list of unique nodes from inodes1 and inodes2,
+! inodes1_u = unique(isort1)
+! inodes2_u = unique(isort2)
+! 2. sort the nodes by coordinates. Now both faces correspond.
+! [coord,k1] = sort(coords(inodes1_u))
+! k1 = inodes1_u(k1)
+! [coord,k2] = sort(coords(inodes2_u))
+! k2 = inodes2_u(k2)
+! 3. set the coordinates on both sides equal to their average
+! coords(k1) = 0.5*( coords(k1)+coords(k2) )
+! coords(k2) = coords(k1)
+ subroutine close_fault_single(fb,nodes_coords,nnodes)
+
+ integer ,intent(in) :: nnodes
+ type(fault_type),intent(in) :: fb
+ double precision,dimension(3,nnodes), intent(inout) :: nodes_coords
+
+ double precision, dimension(3) :: xyz_1, xyz_2
+ double precision, dimension(3) :: xyz
+
+ double precision :: dist
+ integer :: iglob1, iglob2, i, j, k1, k2
+ logical :: found_it
+
+ do e = 1,fb%nspec
+ do k2=1,4
+ iglob2 = fb%inodes2(k2,e)
+ found_it = .false.
+ xyz_2 = nodes_coords(:,iglob2)
+
+ do j = 1,fb%nspec
+ do k1=1,4
+ iglob1 = fb%inodes1(k1,e)
+ xyz_1 = nodes_coords(:,iglob1)
+ xyz = xyz_2-xyz_1
+ dist = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
+ dist = sqrt(dist)
+
+ !jpa: Closing nodes that are already closed is not a problem
+ !jpa: I process them again to leave the loop as early as possible
+ !jpa: and to test if a facing node was found (See below).
+
+ if (dist <= FAULT_GAP_TOLERANCE) then
+ xyz = (xyz_1 + xyz_2)*0.5d0
+ nodes_coords(:,iglob2) = xyz
+ nodes_coords(:,iglob1) = xyz
+ found_it = .true.
+ exit
+ endif
+
+ enddo
+ if (found_it) exit
+ enddo
+
+ if (.not.found_it) then
+ ! If the fault is very complicated (non-planar) the meshes of the two fault sides
+ ! can be very unstructured and might not match (unless you enforce it in CUBIT).
+ ! That is unlikely because the fault gap is tiny, but we never know.
+ stop 'Inconsistent fault mesh: corresponding node in the other fault face was not found'
+ endif
+
+ enddo
+ enddo
+
+ end subroutine close_fault_single
+
+! ---------------------------------------------------------------------------------------------------
+ !--------------------------------------------------
+ ! Repartitioning : two coupled faultside1/side2 elements are transfered to the same partition
+ !--------------------------------------------------
+
+ Subroutine fault_collect_elements(nelmnts,nnodes,elmnts, &
+ sup_neighbour,esize,nsize,nproc,part)
+
+! INPUTS
+ integer(long), intent(in) :: nelmnts,nsize,sup_neighbour
+ integer, dimension(0:esize*nelmnts-1),intent(in) :: elmnts
+ integer, intent(in) :: nnodes, nproc, esize
+! OUTPUTS :
+ integer, dimension(0:nelmnts-1),intent(inout) :: part
+! VARIABLES:
+ logical, dimension(nelmnts) :: is_on_fault
+ integer :: iflt
+
+ is_on_fault = .false.
+ do iflt=1,size(faults)
+ is_on_fault(faults(iflt)%ispec1) = .true.
+ is_on_fault(faults(iflt)%ispec2) = .true.
+ end do
+ call fault_repartition (nelmnts, nnodes, elmnts, sup_neighbour, nsize, &
+ nproc, part, is_on_fault,esize)
+
+ end Subroutine fault_collect_elements
+
+! ---------------------------------------------------------------------------------------------------
+
+ Subroutine fault_repartition (nelmnts, nnodes, elmnts, sup_neighbour, nsize, &
+ nproc, part, is_on_fault, esize)
+
+! INDIVIDUAL FAULT REPARTITION
+
+! part : iproc number of processor partioned. It will altered patching fault elements into the same partion.
+! Part, once is altered , will be input for write_partition_database.
+
+!INPUTS
+ integer(long), intent(in) :: nelmnts,sup_neighbour,nsize
+ integer, intent(in) :: nnodes, nproc, esize
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+ logical , dimension(nelmnts), intent(in) :: is_on_fault
+!OUTPUTS :
+ integer, dimension(0:nelmnts-1), intent(inout) :: part
+
+!LOCAL VARIABLES :
+ integer, dimension(0:nelmnts) :: xadj
+ integer, dimension(0:sup_neighbour*nelmnts-1) :: adjncy
+ integer, dimension(0:nnodes-1) :: nnodes_elmnts
+ integer, dimension(0:nsize*nnodes-1) :: nodes_elmnts
+ integer :: max_neighbour
+
+!SHILDING
+ integer :: i,j, ipart,nproc_null,nproc_null_final
+ integer :: el, el_1, el_2, k1, k2
+ logical :: is_repartitioned
+ integer, dimension(:), allocatable :: elem_proc_null
+
+ ! downloading processor 0
+ nproc_null = count( part == 0 )
+
+ print*, 'Elements proc = 0 redistributed in [{nproc}- nproc0] :'
+ print*, nproc_null
+
+ if ( nproc_null /= 0 ) then
+
+ allocate(elem_proc_null(nproc_null))
+ ! Filling up proc = 0 elements
+ nproc_null = 0
+ do i = 1,nelmnts
+ if ( part(i) == 0 ) then
+ nproc_null = nproc_null + 1
+ elem_proc_null(nproc_null) = i
+ end if
+ end do
+ ! Redistributing proc-0 elements on the rest of processors
+ !jpa: why do this? does it always help balancing ?
+ !pgb: Yes, bulk elements in processor 0 are taken out and redistributed.
+ !pgb: leaving more space for fault elements.
+ !jpa: But if the number of fault elements is much smaller than nproc_null
+ ! we will end up with a very UNbalanced proc 0 !
+ ipart=0
+ do i = 1, nproc_null
+ if ( ipart == nproc ) ipart = 0
+ ipart = ipart +1
+ part(elem_proc_null(i)) = ipart
+ end do
+
+ endif
+ call mesh2dual_ncommonnodes_fault(nelmnts, nnodes, nsize, sup_neighbour, &
+ elmnts, xadj, adjncy, nnodes_elmnts, &
+ nodes_elmnts, max_neighbour, 4, esize)
+
+ ! coupled elements
+ ! ---------------
+ ! Allocating neighbours with shared fault faces.
+ ! This is done in order to not break MPI interfaces.
+ ! Dirty implementation.
+ ! Fault faces are shield by double coarse neighbours.
+ !
+ ! 1 2
+ ! 1 1 2 2
+ ! 1 1 [1] [2] 2 2
+ ! 1 1 2 2
+ ! 1 2
+
+ ! Allocating elements with double shield layer
+ print *, "Fault shield double-layer :"
+ do el = 0, nelmnts-1
+ if ( is_on_fault(el+1) ) then
+ part(el) = 0
+ do k1 = xadj(el), xadj(el+1) - 1
+ el_1 = adjncy(k1)
+ part(el_1) = 0
+ do k2 = xadj(el_1), xadj(el_1+1) - 1
+ el_2 = adjncy(k2)
+ part(el_2) = 0
+ enddo
+ enddo
+ endif
+ enddo
+
+ nproc_null_final = count( part == 0 )
+ print *, nproc_null_final
+
+ end subroutine fault_repartition
+
+! ---------------------------------------------------------------------------------------------------
+! See subroutine write_boundaries_database in part_decompose_mesh_SCOTCH.f90
+!
+! File format:
+! one block for each fault
+! first line of each block = number of fault elements in this processor
+! next lines: #id_(element containing the face) #id_node1_face .. #id_node4_face
+! first for all faces on side 1, then side 2
+
+ subroutine write_fault_database(IIN_database, iproc, nelmnts, glob2loc_elmnts, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part)
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: iproc
+ integer(long), intent(in) :: nelmnts
+ integer, dimension(0:nelmnts-1), intent(in) :: part
+ integer, dimension(0:nelmnts-1), intent(in) :: glob2loc_elmnts
+ integer, dimension(0:), intent(in) :: glob2loc_nodes_nparts
+ integer, dimension(0:), intent(in) :: glob2loc_nodes_parts
+ integer, dimension(0:), intent(in) :: glob2loc_nodes
+
+ integer :: i,j,k,iflt,e
+ integer :: nspec_fault_1,nspec_fault_2
+ integer :: loc_nodes(4),inodes(4)
+
+ do iflt=1,size(faults)
+
+ ! get number of fault elements in this partition
+ nspec_fault_1 = count( part(faults(iflt)%ispec1-1) == iproc )
+ nspec_fault_2 = count( part(faults(iflt)%ispec2-1) == iproc )
+ if (nspec_fault_1 /= nspec_fault_2) then
+ print *, 'Fault # ',iflt,', proc # ',iproc
+ print *, ' ispec1 : ', nspec_fault_1
+ print *, ' ispec2 : ', nspec_fault_2
+ print *, 'Fatal error: Number of fault elements do not coincide. Abort.'
+ stop
+ end if
+ write(IIN_database,*) nspec_fault_1
+
+ ! if no fault element in this partition, move to next fault
+ if (nspec_fault_1==0) cycle
+
+ ! export fault element data, side 1
+ do i=1,faults(iflt)%nspec
+ e = faults(iflt)%ispec1(i)
+ if (part(e-1) == iproc) then
+ inodes = faults(iflt)%inodes1(:,i)
+ do k=1,4
+ do j = glob2loc_nodes_nparts(inodes(k)-1), glob2loc_nodes_nparts(inodes(k))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_nodes(k) = glob2loc_nodes(j) + 1
+ end if
+ end do
+ end do
+ write(IIN_database,*) glob2loc_elmnts(e-1)+1, loc_nodes
+ end if
+ enddo
+
+ ! export fault element data, side 2
+ do i=1,faults(iflt)%nspec
+ e = faults(iflt)%ispec2(i)
+ if(part(e-1) == iproc) then
+ inodes = faults(iflt)%inodes2(:,i)
+ do k=1,4
+ do j = glob2loc_nodes_nparts(inodes(k)-1), glob2loc_nodes_nparts(inodes(k))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_nodes(k) = glob2loc_nodes(j)+1
+ end if
+ end do
+ end do
+ write(IIN_database,*) glob2loc_elmnts(e-1)+1, loc_nodes
+ end if
+ enddo
+
+ enddo
+
+ end subroutine write_fault_database
+
+
+! ---------------------------------------------------------------------------------------------------
+! Creating dual graph (adjacency is defined by 'ncommonnodes' between two elements).
+! Taken from part_decomposition_SCOTCH.f90 routine.
+!----------------------------------------------------------------------------------------------------
+ subroutine mesh2dual_ncommonnodes_fault(nelmnts, nnodes, nsize, sup_neighbour, elmnts,&
+ xadj, adjncy, &
+ nnodes_elmnts, nodes_elmnts, &
+ max_neighbour, ncommonnodes,esize)
+
+ integer(long), intent(in) :: nelmnts
+ integer, intent(in) :: nnodes
+ integer(long), intent(in) :: nsize
+ integer(long), intent(in) :: sup_neighbour
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+
+ integer, dimension(0:nelmnts) :: xadj
+ integer, dimension(0:sup_neighbour*nelmnts-1) :: adjncy
+ integer, dimension(0:nnodes-1) :: nnodes_elmnts
+ integer, dimension(0:nsize*nnodes-1) :: nodes_elmnts
+ integer, intent(out) :: max_neighbour
+ integer, intent(in) :: ncommonnodes,esize
+
+ ! local parameters
+ integer :: i, j, k, l, m, nb_edges
+ logical :: is_neighbour
+ integer :: num_node, n
+ integer :: elem_base, elem_target
+ integer :: connectivity
+
+
+ ! initializes
+ xadj(:) = 0
+ adjncy(:) = 0
+ nnodes_elmnts(:) = 0
+ nodes_elmnts(:) = 0
+ nb_edges = 0
+
+ ! list of elements per node
+ do i = 0, esize*nelmnts-1
+ nodes_elmnts(elmnts(i)*nsize+nnodes_elmnts(elmnts(i))) = i/esize
+ nnodes_elmnts(elmnts(i)) = nnodes_elmnts(elmnts(i)) + 1
+ end do
+
+ ! checking which elements are neighbours ('ncommonnodes' criteria)
+ do j = 0, nnodes-1
+ do k = 0, nnodes_elmnts(j)-1
+ do l = k+1, nnodes_elmnts(j)-1
+
+ connectivity = 0
+ elem_base = nodes_elmnts(k+j*nsize)
+ elem_target = nodes_elmnts(l+j*nsize)
+ do n = 1, esize
+ num_node = elmnts(esize*elem_base+n-1)
+ do m = 0, nnodes_elmnts(num_node)-1
+ if ( nodes_elmnts(m+num_node*nsize) == elem_target ) then
+ connectivity = connectivity + 1
+ end if
+ end do
+ end do
+
+ if ( connectivity >= ncommonnodes) then
+
+ is_neighbour = .false.
+
+ do m = 0, xadj(nodes_elmnts(k+j*nsize))
+ if ( .not.is_neighbour ) then
+ if ( adjncy(nodes_elmnts(k+j*nsize)*sup_neighbour+m) == nodes_elmnts(l+j*nsize) ) then
+ is_neighbour = .true.
+
+ end if
+ end if
+ end do
+ if ( .not.is_neighbour ) then
+ adjncy(nodes_elmnts(k+j*nsize)*sup_neighbour+xadj(nodes_elmnts(k+j*nsize))) = nodes_elmnts(l+j*nsize)
+
+ xadj(nodes_elmnts(k+j*nsize)) = xadj(nodes_elmnts(k+j*nsize)) + 1
+ if (xadj(nodes_elmnts(k+j*nsize))>sup_neighbour) stop 'ERROR : too much neighbours per element, modify the mesh.'
+
+ adjncy(nodes_elmnts(l+j*nsize)*sup_neighbour+xadj(nodes_elmnts(l+j*nsize))) = nodes_elmnts(k+j*nsize)
+
+ xadj(nodes_elmnts(l+j*nsize)) = xadj(nodes_elmnts(l+j*nsize)) + 1
+ if (xadj(nodes_elmnts(l+j*nsize))>sup_neighbour) stop 'ERROR : too much neighbours per element, modify the mesh.'
+ end if
+ end if
+ end do
+ end do
+ end do
+
+ max_neighbour = maxval(xadj)
+
+ ! making adjacency arrays compact (to be used for partitioning)
+ do i = 0, nelmnts-1
+ k = xadj(i)
+ xadj(i) = nb_edges
+ do j = 0, k-1
+ adjncy(nb_edges) = adjncy(i*sup_neighbour+j)
+ nb_edges = nb_edges + 1
+ end do
+ end do
+
+ xadj(nelmnts) = nb_edges
+
+
+ end subroutine mesh2dual_ncommonnodes_fault
+
+!-----------
+! subroutine write_fault_database_iface(IIN_database, iproc, nelmnts, &
+! glob2loc_elmnts, part)
+!
+! integer, intent(in) :: IIN_database
+! integer, intent(in) :: iproc
+! integer(long), intent(in) :: nelmnts
+! integer, dimension(0:nelmnts-1), intent(in) :: part,glob2loc_elmnts
+!
+! integer, dimension(:), allocatable :: ispec1, ispec2, iface1, iface2
+! integer :: i,iflt,ispec_fault
+! integer :: nspec_fault_1,nspec_fault_2,nspec_fault
+! integer :: k1,k2
+!
+! do iflt=1,size(faults)
+!
+! ! check number of fault elements in this partition
+! nspec_fault_1 = count( part(faults(iflt)%ispec1-1) == iproc )
+! nspec_fault_2 = count( part(faults(iflt)%ispec2-1) == iproc )
+! print *, 'Fault # ',iflt
+! print *, ' ispec1 : ', nspec_fault_1
+! print *, ' ispec2 : ', nspec_fault_2
+! if (nspec_fault_1 /= nspec_fault_2) then
+! print *, 'Fatal error: Number of fault elements on ',iproc,' do not coincide. Abort.'
+! stop
+! end if
+! nspec_fault = nspec_fault_1
+!
+! write(IIN_database,*) nspec_fault
+!
+! if (nspec_fault==0) cycle
+!
+! allocate(ispec1(nspec_fault))
+! allocate(iface1(nspec_fault))
+! allocate(ispec2(nspec_fault))
+! allocate(iface2(nspec_fault))
+! k1 = 0
+! k2 = 0
+! do i = 1,faults(iflt)%nspec
+! if ( part(faults(iflt)%ispec1(i)-1) == iproc ) then
+! k1 = k1 + 1
+! ispec1(k1)=glob2loc_elmnts(faults(iflt)%ispec1(i))
+! iface1(k1)=faults(iflt)%iface1(i)
+! endif
+! if ( part(faults(iflt)%ispec2(i)-1) == iproc ) then
+! k2 = k2 + 1
+! ispec2(k2)=glob2loc_elmnts(faults(iflt)%ispec2(i))
+! iface2(k2)=faults(iflt)%iface2(i)
+! endif
+! enddo
+!
+! ! Writes ispec1 , ispec2 , iface1 , iface2
+! do i = 1,nspec_fault
+! write(IIN_database,*) ispec1(i), ispec2(i), iface1(i), iface2(i)
+! enddo
+! ! NOTE: the solver does not need ispec1 and ispec2 to be facing each other across the fault
+! deallocate(ispec1,ispec2,iface1,iface2)
+!
+! enddo
+!
+! end subroutine write_fault_database_iface
+
+
+
+end module fault_scotch
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2011-05-16 00:11:00 UTC (rev 18378)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2011-05-16 02:06:14 UTC (rev 18379)
@@ -189,20 +189,19 @@
nodes_coords(:,iglob1) = xyz
found_it = .true.
exit
- endif
+ endif
+
enddo
+ if (found_it) exit
+ enddo
- if (found_it) then
- exit
- else
+ enddo
+ ! if (.not.found_it) then
! If the fault is very complicated (non-planar) the meshes of the two fault sides
! can be very unstructured and might not match (unless you enforce it in CUBIT).
! That is unlikely because the fault gap is tiny, but we never know.
- stop 'Inconsistent fault mesh: corresponding node in the other fault face was not found'
- endif
- enddo
-
- enddo
+ ! stop 'Inconsistent fault mesh: corresponding node in the other fault face was not found'
+ ! endif
enddo
end subroutine close_fault_single
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/lap_to_brutus
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/lap_to_brutus (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/lap_to_brutus 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,6 @@
+rsync -avuze ssh --delete . gpercy at brutus.ethz.ch:~/FAULT_SOURCE/
+
+## v : verbose , u : preserve newer existing files from overwritting.
+## b : backup existing files before overwritting.
+## e : allows to do remote login ssh.
+## delete : delete different files
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/lap_to_brutus
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/bigstar/run.xgenerate_databases
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/bigstar/run.xgenerate_databases (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/bigstar/run.xgenerate_databases 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+#rm DATABASES_MPI/*.*
+
+# d=`date`
+# echo "Starting compilation $d"
+# make generate_databases
+# d=`date`
+# echo "Finished compilation $d"
+
+
+# compute total number of nodes needed
+NPROC=`grep NPROC DATA/Par_file | cut -d = -f 2`
+
+# total number of nodes is the product of the values read
+
+echo "Submitting job"
+echo $NPROC
+echo "running mesh"
+LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/lib mpirun -np $NPROC ./xgenerate_databases
+echo "database done"
\ No newline at end of file
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/bigstar/run.xgenerate_databases
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/bigstar/run.xspecfem3d
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/bigstar/run.xspecfem3d (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/bigstar/run.xspecfem3d 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,20 @@
+#!/bin/bash
+
+#rm OUTPUT_FILES/*
+
+# d=`date`
+# echo "Starting compilation $d"
+# make specfem3D
+# d=`date`
+# echo "Finished compilation $d"
+
+# compute total number of nodes needed
+NPROC=`grep NPROC DATA/Par_file | cut -d = -f 2`
+
+# total number of nodes is the product of the values read
+
+echo "Submitting job"
+echo $NPROC
+echo "running solver"
+LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/lib mpirun -np $NPROC ./xspecfem3D &
+echo "solver done"
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/bigstar/run.xspecfem3d
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/go_generate_databases_lsf.bash
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/go_generate_databases_lsf.bash (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/go_generate_databases_lsf.bash 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,35 @@
+#!/bin/bash
+
+## job name and output file
+#BSUB -J go_generate_databases
+#BSUB -o OUTPUT_FILES/%J.o
+#BSUB -e OUTPUT_FILES/%J.e
+
+###########################################################
+
+if [ -z $USER ]; then
+ echo "could not run go_mesher_...bash as no USER env is set"
+ exit 2
+fi
+
+# script to run the mesher and the solver
+# read DATA/Par_file to get information about the run
+# compute total number of nodes needed
+NPROC=`grep NPROC DATA/Par_file | cut -d = -f 2 `
+
+# total number of nodes is the product of the values read
+numnodes=$NPROC
+
+cp DATA/Par_file OUTPUT_FILES/
+
+# obtain lsf job information
+cat $LSB_DJOB_HOSTFILE > OUTPUT_FILES/compute_nodes
+echo "$LSB_JOBID" > OUTPUT_FILES/jobid
+
+echo starting MPI mesher on $numnodes processors
+echo " "
+
+sleep 2
+mpirun $PWD/xgenerate_databases
+
+echo "done "
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/go_generate_databases_lsf.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/go_specfem3D_lsf.bash
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/go_specfem3D_lsf.bash (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/go_specfem3D_lsf.bash 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,34 @@
+#!/bin/bash
+
+## job name and output file
+#BSUB -J go_generate_databases
+#BSUB -o OUTPUT_FILES/%J.o
+#BSUB -e OUTPUT_FILES/%J.e
+
+###########################################################
+
+if [ -z $USER ]; then
+ echo "could not run go_mesher_...bash as no USER env is set"
+ exit 2
+fi
+
+# script to run the mesher and the solver
+# read DATA/Par_file to get information about the run
+# compute total number of nodes needed
+NPROC=`grep NPROC DATA/Par_file | cut -d = -f 2 `
+
+# total number of nodes is the product of the values read
+numnodes=$NPROC
+
+cp DATA/Par_file OUTPUT_FILES/
+
+# obtain lsf job information
+cat $LSB_DJOB_HOSTFILE > OUTPUT_FILES/compute_nodes
+echo "$LSB_JOBID" > OUTPUT_FILES/jobid
+
+echo starting MPI mesher on $numnodes processors
+echo " "
+
+sleep 2
+mpirun $PWD/xspecfem3D
+echo "done "
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/go_specfem3D_lsf.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/run_generate_database.bash
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/run_generate_database.bash (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/run_generate_database.bash 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,34 @@
+#!/bin/bash
+
+# this script launches the specfem simulation
+# Percy Galvez, ETH-zurich.Brutus.. 2010.
+
+# use the normal queue unless otherwise directed
+
+queue=""
+if [ $# -eq 1 ]; then
+ echo "Setting the queue to $1"
+ queue="-q $1"
+fi
+
+d=`date`
+echo "Starting compilation $d"
+
+make clean
+make
+
+d=`date`
+
+echo "Finished compilation $d"
+
+# total number of nodes is the product of the values read
+numnodes=`grep NPROC DATA/Par_file | cut -d = -f 2`
+
+echo "Submitting job"
+echo "bsub $queue -R ib -n $numnodes -W 60 < go_generate_databases_lsf.bash"
+
+# Intel version of open MPI needs infiniBand libraries (librdmacm.so),even
+#therefore nodes with this characteristic are required, in that sense option
+# -R ib will take infiniband nodes.
+
+bsub -R ib -n $numnodes -W 60 < go_generate_databases_lsf.bash
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/run_generate_database.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/run_specfem3d.bash
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/run_specfem3d.bash (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/run_specfem3d.bash 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,33 @@
+#!/bin/bash
+
+# this script launches the specfem simulation
+# Percy Galvez, ETH-zurich.Brutus.. 2010.
+
+# use the normal queue unless otherwise directed
+
+queue=""
+if [ $# -eq 1 ]; then
+ echo "Setting the queue to $1"
+ queue="-q $1"
+fi
+
+d=`date`
+echo "Starting compilation $d"
+make clean
+make
+
+d=`date`
+
+echo "Finished compilation $d"
+
+# total number of nodes is the product of the values read
+numnodes=`grep NPROC DATA/Par_file | cut -d = -f 2`
+
+echo "Submitting job"
+echo "bsub $queue -R ib -n $numnodes -W 360 < go_specfem3D_lsf.bash"
+
+# Intel version of open MPI needs infiniBand libraries (librdmacm.so),even
+#therefore nodes with this characteristic are required, in that sense option
+# -R ib will take infiniband nodes.
+
+bsub -R ib -n $numnodes -W 360 < go_specfem3D_lsf.bash
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/brutus/run_specfem3d.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/go_generate_databases_lsf.bash
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/go_generate_databases_lsf.bash (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/go_generate_databases_lsf.bash 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,35 @@
+#!/bin/bash
+
+## job name and output file
+#BSUB -J go_generate_databases
+#BSUB -o OUTPUT_FILES/%J.o
+#BSUB -e OUTPUT_FILES/%J.e
+
+###########################################################
+
+if [ -z $USER ]; then
+ echo "could not run go_mesher_...bash as no USER env is set"
+ exit 2
+fi
+
+# script to run the mesher and the solver
+# read DATA/Par_file to get information about the run
+# compute total number of nodes needed
+NPROC=`grep NPROC DATA/Par_file | cut -d = -f 2 `
+
+# total number of nodes is the product of the values read
+numnodes=$NPROC
+
+cp DATA/Par_file OUTPUT_FILES/
+
+# obtain lsf job information
+cat $LSB_DJOB_HOSTFILE > OUTPUT_FILES/compute_nodes
+echo "$LSB_JOBID" > OUTPUT_FILES/jobid
+
+echo starting MPI mesher on $numnodes processors
+echo " "
+
+sleep 2
+mpirun $PWD/xgenerate_databases
+
+echo "done "
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/go_generate_databases_lsf.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/go_specfem3D_lsf.bash
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/go_specfem3D_lsf.bash (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/go_specfem3D_lsf.bash 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,34 @@
+#!/bin/bash
+
+## job name and output file
+#BSUB -J go_generate_databases
+#BSUB -o OUTPUT_FILES/%J.o
+#BSUB -e OUTPUT_FILES/%J.e
+
+###########################################################
+
+if [ -z $USER ]; then
+ echo "could not run go_mesher_...bash as no USER env is set"
+ exit 2
+fi
+
+# script to run the mesher and the solver
+# read DATA/Par_file to get information about the run
+# compute total number of nodes needed
+NPROC=`grep NPROC DATA/Par_file | cut -d = -f 2 `
+
+# total number of nodes is the product of the values read
+numnodes=$NPROC
+
+cp DATA/Par_file OUTPUT_FILES/
+
+# obtain lsf job information
+cat $LSB_DJOB_HOSTFILE > OUTPUT_FILES/compute_nodes
+echo "$LSB_JOBID" > OUTPUT_FILES/jobid
+
+echo starting MPI mesher on $numnodes processors
+echo " "
+
+sleep 2
+mpirun $PWD/xspecfem3D
+echo "done "
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/go_specfem3D_lsf.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_generate_database.bash
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_generate_database.bash (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_generate_database.bash 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,34 @@
+#!/bin/bash
+
+# this script launches the specfem simulation
+# Percy Galvez, ETH-zurich.Brutus.. 2010.
+
+# use the normal queue unless otherwise directed
+
+queue=""
+if [ $# -eq 1 ]; then
+ echo "Setting the queue to $1"
+ queue="-q $1"
+fi
+
+d=`date`
+echo "Starting compilation $d"
+
+make clean
+make
+
+d=`date`
+
+echo "Finished compilation $d"
+
+# total number of nodes is the product of the values read
+numnodes=`grep NPROC DATA/Par_file | cut -d = -f 2`
+
+echo "Submitting job"
+echo "bsub $queue -R ib -n $numnodes -W 60 < go_generate_databases_lsf.bash"
+
+# Intel version of open MPI needs infiniBand libraries (librdmacm.so),even
+#therefore nodes with this characteristic are required, in that sense option
+# -R ib will take infiniband nodes.
+
+bsub -R ib -n $numnodes -W 60 < go_generate_databases_lsf.bash
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_generate_database.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_lsf.bash
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_lsf.bash (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_lsf.bash 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,33 @@
+#!/bin/bash
+
+# this script launches the specfem simulation
+# Percy Galvez, ETH-zurich.Brutus.. 2010.
+
+# use the normal queue unless otherwise directed
+
+queue=""
+if [ $# -eq 1 ]; then
+ echo "Setting the queue to $1"
+ queue="-q $1"
+fi
+
+d=`date`
+echo "Starting compilation $d"
+make clean
+make
+
+d=`date`
+
+echo "Finished compilation $d"
+
+# total number of nodes is the product of the values read
+numnodes=`grep NPROC DATA/Par_file | cut -d = -f 2`
+
+echo "Submitting job"
+echo "bsub $queue -R ib -n $numnodes -W 360 < go_specfem3D_lsf.bash"
+
+# Intel version of open MPI needs infiniBand libraries (librdmacm.so),even
+#therefore nodes with this characteristic are required, in that sense option
+# -R ib will take infiniband nodes.
+
+bsub -R ib -n $numnodes -W 360 < go_specfem3D_lsf.bash
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_lsf.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_specfem3d.bash
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_specfem3d.bash (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_specfem3d.bash 2011-05-16 02:06:14 UTC (rev 18379)
@@ -0,0 +1,33 @@
+#!/bin/bash
+
+# this script launches the specfem simulation
+# Percy Galvez, ETH-zurich.Brutus.. 2010.
+
+# use the normal queue unless otherwise directed
+
+queue=""
+if [ $# -eq 1 ]; then
+ echo "Setting the queue to $1"
+ queue="-q $1"
+fi
+
+d=`date`
+echo "Starting compilation $d"
+make clean
+make
+
+d=`date`
+
+echo "Finished compilation $d"
+
+# total number of nodes is the product of the values read
+numnodes=`grep NPROC DATA/Par_file | cut -d = -f 2`
+
+echo "Submitting job"
+echo "bsub $queue -R ib -n $numnodes -W 360 < go_specfem3D_lsf.bash"
+
+# Intel version of open MPI needs infiniBand libraries (librdmacm.so),even
+#therefore nodes with this characteristic are required, in that sense option
+# -R ib will take infiniband nodes.
+
+bsub -R ib -n $numnodes -W 360 < go_specfem3D_lsf.bash
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/run/run_specfem3d.bash
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2011-05-16 00:11:00 UTC (rev 18378)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2011-05-16 02:06:14 UTC (rev 18379)
@@ -72,7 +72,7 @@
iface3_corner_ijk,iface4_corner_ijk, &
iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
- public :: fault_read_input, fault_setup, fault_save_arrays_test, fault_save_arrays, fault_db, &
+ public :: fault_read_input, fault_setup, fault_save_arrays_test, fault_save_arrays, &
nodes_coords_open, nnodes_coords_open, ANY_FAULT_IN_THIS_PROC
contains
@@ -191,7 +191,7 @@
call setup_Kelvin_Voigt_eta(fault_db(iflt),nspec)
- call save_fault_xyzcoord_ibulk(fault_db(iflt),ibool,nspec)
+ call save_fault_xyzcoord_ibulk(fault_db(iflt))
call setup_normal_jacobian(fault_db(iflt),ibool,nspec,nglob,myrank)
@@ -201,61 +201,6 @@
!=============================================================================================================
-! We only close *store_dummy. Closing *store is more difficult.
-! Fortunately only *store_dummy is needed to compute jacobians and normals
-subroutine close_fault(fdb) !,xstore,ystore,zstore,nspec)
-
- use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy
-
- type(fault_db_type), intent(inout) :: fdb
-! integer, intent(in) :: nspec
-! double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(inout) :: xstore,ystore,zstore
-
- integer :: i,K1,K2
-
- do i=1,fdb%nglob
- K1 = fdb%ibulk1(i)
- K2 = fdb%ibulk2(i)
- xstore_dummy(K1) = 0.5d0*( xstore_dummy(K1) + xstore_dummy(K2) )
- xstore_dummy(K2) = xstore_dummy(K1)
- ystore_dummy(K1) = 0.5d0*( ystore_dummy(K1) + ystore_dummy(K2) )
- ystore_dummy(K2) = ystore_dummy(K1)
- zstore_dummy(K1) = 0.5d0*( zstore_dummy(K1) + zstore_dummy(K2) )
- zstore_dummy(K2) = zstore_dummy(K1)
- enddo
-
-! Closing *store is not easy ...
-! The following would work only if ispec1 and ispec2 were always facing each other
-! and if igll corresponded to the same node in both element faces.
-! That would require reordering isepcs, e.g. sorting them by coordinate of the middle point of the face,
-! and then, within each face, searching the nearest node (by coordinate)
-!
-! do e=1,fdb%nspec
-! e1 = fdb%ispec1(e)
-! e2 = fdb%ispec1(e)
-! igll = 0
-! do i=1,NGLLX
-! do j=1,NGLLX
-! igll = igll + 1
-! i1=fdb%ijk1(1,igll,e)
-! j1=fdb%ijk1(2,igll,e)
-! k1=fdb%ijk1(3,igll,e)
-! i2=fdb%ijk2(1,igll,e)
-! j2=fdb%ijk2(2,igll,e)
-! k2=fdb%ijk2(3,igll,e)
-! xstore(i1,j1,k1,e1) = 0.5d0*( xstore(i1,j1,k1,e1) + xstore(i2,j2,k2,e2) )
-! xstore(i2,j2,k2,e2) = xstore(i1,j1,k1,e1)
-! ystore(i1,j1,k1,e1) = 0.5d0*( ystore(i1,j1,k1,e1) + ystore(i2,j2,k2,e2) )
-! ystore(i2,j2,k2,e2) = ystore(i1,j1,k1,e1)
-! zstore(i1,j1,k1,e1) = 0.5d0*( zstore(i1,j1,k1,e1) + zstore(i2,j2,k2,e2) )
-! zstore(i2,j2,k2,e2) = zstore(i1,j1,k1,e1)
-! enddo
-! enddo
-! enddo
-
-subroutine close_fault
-
-!=============================================================================================================
! looks for i,j,k indices of GLL points on boundary face
! determines element face by given CUBIT corners
! sets face id of reference element associated with this face
@@ -303,7 +248,7 @@
type(fault_db_type), intent(inout) :: fdb
- integer :: e,i,j,igll, iface_ref1, iface_ref2
+ integer :: e,i,j,igll
integer :: ijk_face1(3,NGLLX,NGLLY), ijk_face2(3,NGLLX,NGLLY)
allocate(fdb%ijk1(3,NGLLX*NGLLY,fdb%nspec))
@@ -470,12 +415,11 @@
!=================================================================================
- subroutine save_fault_xyzcoord_ibulk(fdb,ibool,nspec)
+ subroutine save_fault_xyzcoord_ibulk(fdb)
use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy
type(fault_db_type), intent(inout) :: fdb
- integer, intent(in) :: nspec, ibool(NGLLX,NGLLY,NGLLZ,nspec)
integer :: K1, K2, i
Modified: seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/tpv5/DATA/FAULT/Par_file_faults.in
===================================================================
--- seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/tpv5/DATA/FAULT/Par_file_faults.in 2011-05-16 00:11:00 UTC (rev 18378)
+++ seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/tpv5/DATA/FAULT/Par_file_faults.in 2011-05-16 02:06:14 UTC (rev 18379)
@@ -1,5 +1,6 @@
1
-1 2 0.00134 ! #tag1 #tag2 #eta(damping)
+0.00134 ! #tag #eta(damping)
+1 ! 1 = dyn 2=kin
1000 ! NTOUT : Number of time steps
100 ! NTSNAP: time interation of snapshots
-1.0e0 ! V_HEALING (-1 : Healing off)
@@ -11,3 +12,4 @@
&DIST2D shape='square', val = 62.0e6, xc = 7500.0e0, yc =0e0, zc= -7500.0e0, r=0e0, l=3000.0e0, lx=0e0, ly=0e0, lz=0e0 /
&SWF mus=10000.0e0,mud=0.525e0,dc=0.4e0,nmus=1,nmud=0,ndc=0 /
&DIST2D shape='rectangle', val = 0.677e0, xc = 0e0, yc =0e0, zc= -7500.0e0, r=0e0, l=0e0, lx=30000.0e0, ly=0.0e0, lz=15000.0e0 /
+&KINPAR kindt=10.0e-3
More information about the CIG-COMMITS
mailing list