[cig-commits] r20783 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: decompose_mesh_SCOTCH src
ampuero at geodynamics.org
ampuero at geodynamics.org
Wed Sep 26 16:15:02 PDT 2012
Author: ampuero
Date: 2012-09-26 16:15:02 -0700 (Wed, 26 Sep 2012)
New Revision: 20783
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/fault_scotch.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
Log:
fixed close_faults, now based on inodes; some clean up
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 2012-09-26 22:01:09 UTC (rev 20782)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-09-26 23:15:02 UTC (rev 20783)
@@ -83,8 +83,6 @@
integer :: q_flag,aniso_flag,idomain_id
double precision :: vp,vs,rho
- logical, parameter :: PARALLEL_FAULT = .false.
-
contains
!----------------------------------------------------------------------------------------------
@@ -391,11 +389,9 @@
if( nspec2D_moho > 0 ) print*, ' nspec2D_moho = ', nspec2D_moho
call read_fault_files(localpath_name)
- if (allocated(faults)) then
+ if (ANY_FAULT) then
call save_nodes_coords(nodes_coords,nnodes)
- call close_faults(nodes_coords,elmnts,nspec,nnodes,esize)
- !SNS
- if (PARALLEL_FAULT) call reorder_fault_elements(nodes_coords,nnodes)
+ call close_faults(nodes_coords,nnodes)
end if
end subroutine read_mesh_files
@@ -533,11 +529,11 @@
! count_def_mat, mat(1,:) , mat_prop, &
! sup_neighbour, nsize, &
! nparts, part, nfaces_coupled, faces_coupled)
- if(.not. PARALLEL_FAULT) then
+ if (PARALLEL_FAULT) then
+ call fault_repartition_parallel (nspec,part,nodes_coords,nnodes)
+ else
! move all fault elements to the same partition (proc=0)
call fault_repartition (nspec, nnodes, elmnts, nsize, nparts, part, esize)
- else
- call fault_repartition_parallel (nspec,part)
endif
! re-partitioning puts moho-surface coupled elements into same partition
@@ -655,8 +651,8 @@
close(15)
- if (.not. allocated(faults)) cycle
- ! write fault database
+ ! write fault database
+ if (ANY_FAULT) then
write(prname, "(i6.6,'_Database_fault')") ipart
open(unit=16,file=outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname,&
status='replace', action='write', form='formatted', iostat = ierr)
@@ -666,31 +662,18 @@
print*,'check if path exists:',outputpath_name(1:len_trim(outputpath_name))
stop
endif
-print*,prname
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)
+ endif
end do
+
print*, 'partitions: '
print*, ' num = ',nparts
print*
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 2012-09-26 22:01:09 UTC (rev 20782)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2012-09-26 23:15:02 UTC (rev 20783)
@@ -12,8 +12,9 @@
type(fault_type), allocatable, save :: faults(:)
double precision, dimension(:,:), allocatable, save :: nodes_coords_open
+ logical, save :: ANY_FAULT = .false.
+ logical, parameter :: PARALLEL_FAULT = .false.
-
integer, parameter :: long = SELECTED_INT_KIND(18)
double precision, parameter :: FAULT_GAP_TOLERANCE = 1.0d0
@@ -21,7 +22,8 @@
! but smaller than the smallest element size
public :: read_fault_files, fault_repartition, fault_repartition_parallel, close_faults, write_fault_database, &
- save_nodes_coords, nodes_coords_open, faults, reorder_fault_elements
+ save_nodes_coords, nodes_coords_open, ANY_FAULT, PARALLEL_FAULT
+
CONTAINS
!==========================================================================================
@@ -39,13 +41,14 @@
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
+ ANY_FAULT = (nbfaults>0)
+ if (.not. ANY_FAULT) return
+ allocate(faults(nbfaults))
+ do iflt = 1 , nbfaults
+ call read_single_fault_file(faults(iflt),iflt,localpath_name)
+ enddo
+
end subroutine read_fault_files
@@ -68,6 +71,13 @@
NTchar(1:len_trim(NTchar))//'.dat'
filename = adjustl(filename)
! reads fault elements and nodes
+ ! File format:
+ ! Line 1:
+ ! number_of_elements_in_side_1 number_of_elements_in_side_2
+ ! Then for all elements that have a face on side 1:
+ ! #id_element #id_global_node1 .. #id_global_node4
+ ! Then the same for side 2.
+ ! Note: element ids start at 1, not 0 (see cubit2specfem3d.py)
open(unit=101, file=filename, status='old', form='formatted', iostat = ier)
if( ier /= 0 ) then
write(6,*) 'Fatal error: file '//filename//' not found'
@@ -82,11 +92,6 @@
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)
enddo
@@ -119,18 +124,15 @@
! ---------------------------------------------------------------------------------------------------
- subroutine close_faults(nodes_coords,elmnts,nelmnts,nnodes,esize)
+ subroutine close_faults(nodes_coords,nnodes)
- integer ,intent(in) :: nnodes, esize
- integer(long), intent(in) :: nelmnts
+ integer, intent(in) :: nnodes
double precision, dimension(3,nnodes), intent(inout) :: nodes_coords
- integer, dimension(esize,nelmnts), intent(in) :: elmnts
- integer :: iflt
+ integer :: i
- do iflt=1,size(faults)
- call close_fault_single(faults(iflt)%ispec1,faults(iflt)%ispec2,&
- elmnts,nodes_coords,nnodes,esize,nelmnts)
+ do i=1,size(faults)
+ call close_fault_single(faults(i),nodes_coords,nnodes)
enddo
end subroutine close_faults
@@ -148,34 +150,30 @@
! 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(ispec1,ispec2,elmnts,nodes_coords,nnodes,esize,nelmnts)
+ subroutine close_fault_single(f,nodes_coords,nnodes)
- integer ,intent(in) :: nnodes, esize
- integer(long), intent(in) :: nelmnts
- integer, dimension(esize,nelmnts), intent(in) :: elmnts
- integer , dimension(:), intent(in) :: ispec1,ispec2
- double precision,dimension(3,nnodes), intent(inout) :: nodes_coords
+ type(fault_type), intent(in) :: f
+ integer, intent(in) :: nnodes
+ double precision, dimension(3,nnodes), intent(inout) :: nodes_coords
- double precision, dimension(3) :: xyz_1, xyz_2
- double precision, dimension(3) :: xyz
+ double precision, dimension(3) :: xyz_1, xyz_2, xyz
double precision :: dist
integer :: iglob1, iglob2, i, j, k1, k2
logical :: found_it
- do i = 1,size(ispec2)
- do k2=1,esize
- iglob2 = elmnts(k2,ispec2(i))
+ do i = 1,f%nspec
+ do k2=1,4
+ iglob2 = f%inodes2(k2,i)
found_it = .false.
xyz_2 = nodes_coords(:,iglob2)
- do j = 1,size(ispec1)
- do k1=1,esize
- iglob1 = elmnts(k1,ispec1(j))
+ do j = 1,f%nspec
+ do k1=1,4
+ iglob1 = f%inodes1(k1,j)
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)
+ dist = sqrt(xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3))
!jpa: Closing nodes that are already closed is not a problem
!jpa: I process them again to leave the loop as early as possible
@@ -202,6 +200,9 @@
end subroutine close_fault_single
!===================================================================================================
+! Lexicographic reordering of fault elements (based on their centroid)
+! to make sure both sides are ordered in the same way
+! and hence elements facing each other have the same index
subroutine reorder_fault_elements(nodes_coords,nnodes)
integer, intent(in) :: nnodes
@@ -430,13 +431,18 @@
end subroutine fault_repartition
! ---------------------------------------------------------------------------------------------------
- subroutine fault_repartition_parallel (nelmnts, part)
+ subroutine fault_repartition_parallel (nelmnts, part, nodes_coords,nnodes)
integer(long), intent(in) :: nelmnts
integer, dimension(0:nelmnts-1), intent(inout) :: part
+ integer, intent(in) :: nnodes
+ double precision,dimension(3,nnodes), intent(in) :: nodes_coords
integer :: i,iflt,e,e1,e2,proc1,proc2
+ ! Reorder both fault sides so that elements facing each other have the same index
+ call reorder_fault_elements(nodes_coords,nnodes)
+
!JPA loop over all faults
!JPA loop over all fault element pairs
!JPA assign both elements to the processor with lowest rank among the pair
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 2012-09-26 22:01:09 UTC (rev 20782)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2012-09-26 23:15:02 UTC (rev 20783)
@@ -55,7 +55,7 @@
integer, save :: nnodes_coords_open
logical, save :: ANY_FAULT_IN_THIS_PROC = .false.
logical, save :: ANY_FAULT = .false.
- logical, save :: PARALLEL_FAULT = .false.
+ logical, parameter :: PARALLEL_FAULT = .false.
! corners indices of reference cube faces
integer,dimension(3,4),parameter :: iface1_corner_ijk = &
More information about the CIG-COMMITS
mailing list