[cig-commits] r18341 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: . CUBIT DATA decompose_mesh_SCOTCH
percygalvez at geodynamics.org
percygalvez at geodynamics.org
Tue May 10 23:12:33 PDT 2011
Author: percygalvez
Date: 2011-05-10 23:12:33 -0700 (Tue, 10 May 2011)
New Revision: 18341
Added:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/DATA/
seismo/3D/FAULT_SOURCE/branches/new_fault_db/DATA/CMTSOLUTION
seismo/3D/FAULT_SOURCE/branches/new_fault_db/DATA/CMTSOLUTION.in
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py
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/decompose_mesh_SCOTCH.f90
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/part_decompose_mesh_SCOTCH.f90
Log:
update 11.05.2011
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-11 03:08:29 UTC (rev 18340)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py 2011-05-11 06:12:33 UTC (rev 18341)
@@ -19,21 +19,21 @@
ybulk=[]
zbulk=[]
-xbulk.append(-21*km) #x1
-xbulk.append(21*km) #x2
-xbulk.append(21*km) #x3
-xbulk.append(-21*km) #x4
+xbulk.append(-18*km) #x1
+xbulk.append(18*km) #x2
+xbulk.append(18*km) #x3
+xbulk.append(-18*km) #x4
-ybulk.append(-16.2*km) #y1
-ybulk.append(-16.2*km) #y2
-ybulk.append(16.2*km) #y3
-ybulk.append(16.2*km) #y4
+ybulk.append(-18*km) #y1
+ybulk.append(-18*km) #y2
+ybulk.append(18*km) #y3
+ybulk.append(18*km) #y4
zbulk=[z_surf]*4
-x.append(-16*km) #x5
+x.append(-9*km) #x5
x.append(0*km) #x6
-x.append(16*km) #x7
+x.append(9*km) #x7
x.append(0*km) #x8
y.append(0.0) #y5
@@ -80,8 +80,12 @@
cubit.cmd("sweep curve 6 vector 0 0 -1 distance "+str(21*km))
#####################################################
+elementsize = 9000
cubit.cmd("imprint all")
cubit.cmd("merge all")
+cubit.cmd("surface 1 size "+str(elementsize))
+cubit.cmd("volume 1 size "+str(elementsize))
+
cubit.cmd("mesh surface 1")
cubit.cmd("mesh volume 1")
#cubit.cmd("unmerge surface 2 3")
@@ -124,7 +128,7 @@
for f in faces:
if dic_quads_fault_u.has_key(f):
nodes=cubit.get_connectivity('Face',f)
- print 'nodes :',nodes[0],nodes[1],nodes[2],nodes[3]
+# print 'nodes :',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)
@@ -135,7 +139,7 @@
for f in faces:
if dic_quads_fault_d.has_key(f):
nodes=cubit.get_connectivity('Face',f)
- print 'nodes :',nodes[0],nodes[1],nodes[2],nodes[3]
+# print 'nodes :',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)
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/DATA/CMTSOLUTION
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/DATA/CMTSOLUTION (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/DATA/CMTSOLUTION 2011-05-11 06:12:33 UTC (rev 18341)
@@ -0,0 +1,13 @@
+PDE 1999 01 01 00 00 00.00 1000000 1000000 -1000000 1 1 dynamic_rupture
+event name: dynamic_rupture
+time shift: 0.0000
+half duration: 0.0
+latitude: 1000000.0
+longitude: 1000000.0
+depth: -1000000.0
+Mrr: 1.258900e+0
+Mtt: 1.258900e+0
+Mpp: 1.258900e+0
+Mrt: 0.000000e0
+Mrp: 0.000000e0
+Mtp: 0.000000e0
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/DATA/CMTSOLUTION.in
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/DATA/CMTSOLUTION.in (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/DATA/CMTSOLUTION.in 2011-05-11 06:12:33 UTC (rev 18341)
@@ -0,0 +1,13 @@
+PDE 1999 01 01 00 00 00.00 1000000 1000000 -1000000 1 1 dynamic_rupture
+event name: dynamic_rupture
+time shift: 0.0000
+half duration: 0.0
+latitude: 1000000.0
+longitude: 1000000.0
+depth: -1000000.0
+Mrr: 1.258900e+0
+Mtt: 1.258900e+0
+Mpp: 1.258900e+0
+Mrt: 0.000000e0
+Mrp: 0.000000e0
+Mtp: 0.000000e0
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-11 03:08:29 UTC (rev 18340)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/compile_all.bash 2011-05-11 06:12:33 UTC (rev 18341)
@@ -1,9 +1,8 @@
#!/bin/bash
make clean
make
-# ./xdecompose_mesh_SCOTCH n input_dir output_dir
# ./xdecompose_mesh_SCOTCH 4 MESH/OUTPUT_FILES ../DATABASES_MPI/
-npart=`grep nparts constants_decompose_mesh_SCOTCH.h | cut -d = -f 2`
-
+#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 'nproc' ../CUBIT/MESH/ ../DATABASES_MPI/
+./xdecompose_mesh_SCOTCH $npart ../CUBIT/MESH/ ../DATABASES_MPI/
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-05-11 03:08:29 UTC (rev 18340)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-05-11 06:12:33 UTC (rev 18341)
@@ -388,9 +388,10 @@
close(98)
if( nspec2D_moho > 0 ) print*, ' nspec2D_moho = ', nspec2D_moho
- call read_fault_files()
+ call read_fault_files(localpath_name)
call close_faults(nodes_coords,elmnts,nspec,nnodes,esize)
-
+ !TEST
+ stop
end subroutine read_mesh_files
!----------------------------------------------------------------------------------------------
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-11 03:08:29 UTC (rev 18340)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2011-05-11 06:12:33 UTC (rev 18341)
@@ -1,8 +1,8 @@
module fault_scotch
implicit none
- private
-
+
+ private
type fault_type
private
integer :: nspec
@@ -10,10 +10,11 @@
integer, dimension(:,:), pointer :: inodes1, inodes2
end type fault_type
- type(fault_type), allocatable, save :: faults(:)
-
+ type(fault_type), allocatable, save :: faults(:)
integer, parameter :: long = SELECTED_INT_KIND(18)
-
+ integer, parameter :: SIZE_REAL = 4 ! single precision
+ integer, parameter :: SIZE_DOUBLE = 8 ! double precision
+ integer, parameter :: CUSTOM_REAL = SIZE_REAL
double precision, parameter :: FAULT_GAP_TOLERANCE = 1e-2_CUSTOM_REAL
public :: read_fault_files, fault_collect_elements, close_faults, write_fault_database
@@ -21,8 +22,9 @@
CONTAINS
!==========================================================================================
- Subroutine read_fault_files
+ 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)
@@ -37,7 +39,7 @@
if (nbfaults>0) then
allocate(faults(nbfaults))
do iflt = 1 , nbfaults
- call read_single_fault_file(faults(iflt),iflt)
+ call read_single_fault_file(faults(iflt),iflt,localpath_name)
enddo
endif
@@ -46,26 +48,33 @@
!---------------------------------------------------------------------------------------------------
- Subroutine read_single_fault_file(f,ifault)
+ Subroutine read_single_fault_file(f,ifault,localpath_name)
-! INPUTS :
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)
-
- open(101,file='../DATA/FAULT/fault_file_'//NTCHAR//'.dat',&
- status='old',action='read',iostat=ier)
+
+ 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
+ print*,'could not open file:',filename
+ stop 'error file open'
+ endif
if( ier == 0 ) then
read(101,*) nspec_side1,nspec_side2
- if (nspec_side1 =/ nspec_side2) stop 'Number of fault nodes at do not match.'
+ if (nspec_side1 /= nspec_side2) stop 'Number of fault nodes at do not match.'
f%nspec = nspec_side1
- read(101,*) f%nspec
allocate(f%ispec1(f%nspec))
allocate(f%ispec2(f%nspec))
allocate(f%inodes1(4,f%nspec))
@@ -77,9 +86,13 @@
! 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)
+ !TEST
+ print*, f%ispec1(e),f%inodes1(:,e)
enddo
do e=1,f%nspec
read(101,*) f%ispec2(e),f%inodes2(:,e)
+ !TEST
+ print*, f%ispec2(e),f%inodes2(:,e)
enddo
! If we ever figure out how to export "ifaces" from CUBIT:
!allocate(f%iface1(f%nspec))
@@ -88,7 +101,7 @@
! read(101,*) f%ispec1(e),f%ispec2(e),f%iface1(e),f%iface2(e)
!enddo
else
- write(6,*) 'Fatal error: file ../DATA/FAULT/fault_file_'//NTCHAR//'.dat not found'
+ write(6,*) 'Fatal error: file '//filename//' not found'
write(6,*) 'Abort'
stop
endif
@@ -149,19 +162,29 @@
iglob2 = elmnts(k2,ispec2(i))
found_it = .false.
xyz_2 = nodes_coords(:,iglob2)
-
+ ! TEST
+ print*,'iblob2:',iglob2
+ print*,'xzy2 :',xyz_2
do j = 1,size(ispec1)
do k1=1,esize
iglob1 = elmnts(k1,ispec1(j))
xyz_1 = nodes_coords(:,iglob1)
-
+ ! TEST
+ print*,'iblob1:',iglob1
+ print*,'xzy1 :',xyz_1
xyz = xyz_2-xyz_1
dist = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
dist = sqrt(dist)
if (dist <= FAULT_GAP_TOLERANCE) then
xyz = (xyz_1 + xyz_2)*0.5_CUSTOM_REAL
+ ! TEST
+ write(6,*) 'iblob2:',iglob2
+ write(6,*) 'iblob1:',iglob1
+ write(6,*) 'xzy1 :',xyz_1
+ write(6,*) 'xzy2 :',xyz_2
+ write(6,*) 'xzy :',xyz
nodes_coords(:,iglob2) = xyz
nodes_coords(:,iglob1) = xyz
found_it = .true.
@@ -178,7 +201,6 @@
end subroutine close_fault_single
! ---------------------------------------------------------------------------------------------------
-
!--------------------------------------------------
! Repartitioning : two coupled faultside1/side2 elements are transfered to the same partition
!--------------------------------------------------
@@ -267,9 +289,9 @@
part(elem_proc_null(i)) = ipart
end do
- call mesh2dual_ncommonnodes(nelmnts, nnodes, nsize, sup_neighbour, &
+ call mesh2dual_ncommonnodes_fault(nelmnts, nnodes, nsize, sup_neighbour, &
elmnts, xadj, adjncy, nnodes_elmnts, &
- nodes_elmnts, max_neighbour, 4)
+ nodes_elmnts, max_neighbour, 4, esize)
! coupled elements
! ---------------
@@ -378,6 +400,111 @@
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)
@@ -440,4 +567,6 @@
!
! end subroutine write_fault_database_iface
+
+
end module fault_scotch
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2011-05-11 03:08:29 UTC (rev 18340)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2011-05-11 06:12:33 UTC (rev 18341)
@@ -519,7 +519,7 @@
!--------------------------------------------------
subroutine write_glob2loc_nodes_database(IIN_database, iproc, npgeo, &
nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
- glob2loc_nodes, nnodes, num_phase,ispec1,ispec2)
+ glob2loc_nodes, nnodes, num_phase)
integer, intent(in) :: IIN_database
integer, intent(in) :: nnodes, iproc, num_phase
More information about the CIG-COMMITS
mailing list