[cig-commits] r18664 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: . decompose_mesh_SCOTCH/devel src
percygalvez at geodynamics.org
percygalvez at geodynamics.org
Tue Jun 28 07:47:58 PDT 2011
Author: percygalvez
Date: 2011-06-28 07:47:58 -0700 (Tue, 28 Jun 2011)
New Revision: 18664
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/devel/fault_scotch.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_movie_shakemap_AVS_DX_GMT.f90
Log:
New Database finished
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT 2011-06-28 12:00:53 UTC (rev 18663)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT 2011-06-28 14:47:58 UTC (rev 18664)
@@ -235,9 +235,8 @@
3. Rupture time files are named Rupture_time_FAULT-id. Their format is 4 columns:
X, Y, Z and rupture time.
-4. The face FAULT plane reference used in the computatinos is the first face defined by the user.
+4. The face FAULT plane reference used in the computatinos takes fault_down as reference.
-
V. POST-PROCESSING AND VISUALIZATION
-------------------------------------
Modified: 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 2011-06-28 12:00:53 UTC (rev 18663)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/devel/fault_scotch.f90 2011-06-28 14:47:58 UTC (rev 18664)
@@ -16,8 +16,7 @@
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.
+ double precision, parameter :: FAULT_GAP_TOLERANCE = 1.0d0
public :: read_fault_files, fault_collect_elements, close_faults, write_fault_database, &
save_nodes_coords, nodes_coords_open, faults
@@ -121,15 +120,18 @@
! ---------------------------------------------------------------------------------------------------
- subroutine close_faults(nodes_coords,nnodes)
+ subroutine close_faults(nodes_coords,elmnts,nelmnts,nnodes,esize)
- integer ,intent(in) :: nnodes
+ integer ,intent(in) :: nnodes, esize
+ integer(long), intent(in) :: nelmnts
double precision, dimension(3,nnodes), intent(inout) :: nodes_coords
+ integer, dimension(esize,nelmnts), intent(in) :: elmnts
integer :: iflt
do iflt=1,size(faults)
- call close_fault_single(faults(iflt),nodes_coords,nnodes)
+ call close_fault_single(faults(iflt)%ispec1,faults(iflt)%ispec2,&
+ elmnts,nodes_coords,nnodes,esize,nelmnts)
enddo
end subroutine close_faults
@@ -147,10 +149,12 @@
! 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)
+ subroutine close_fault_single(ispec1,ispec2,elmnts,nodes_coords,nnodes,esize,nelmnts)
- integer ,intent(in) :: nnodes
- type(fault_type),intent(in) :: fb
+ 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
double precision, dimension(3) :: xyz_1, xyz_2
@@ -160,15 +164,15 @@
integer :: iglob1, iglob2, i, j, k1, k2
logical :: found_it
- do e = 1,fb%nspec
- do k2=1,4
- iglob2 = fb%inodes2(k2,e)
+ do i = 1,size(ispec2)
+ do k2=1,esize
+ iglob2 = elmnts(k2,ispec2(i))
found_it = .false.
xyz_2 = nodes_coords(:,iglob2)
- do j = 1,fb%nspec
- do k1=1,4
- iglob1 = fb%inodes1(k1,e)
+ do j = 1,size(ispec1)
+ do k1=1,esize
+ iglob1 = elmnts(k1,ispec1(j))
xyz_1 = nodes_coords(:,iglob1)
xyz = xyz_2-xyz_1
dist = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
@@ -190,14 +194,13 @@
if (found_it) exit
enddo
- if (.not.found_it) then
+ 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
+ ! stop 'Inconsistent fault mesh: corresponding node in the other fault face was not found'
+ ! endif
enddo
end subroutine close_fault_single
@@ -257,7 +260,7 @@
!SHILDING
integer :: i,j, ipart,nproc_null,nproc_null_final
- integer :: el, el_1, el_2, k1, k2
+ integer :: el, el_1, el_2, k1, k2, k,e,iflt,inode
logical :: is_repartitioned
integer, dimension(:), allocatable :: elem_proc_null
@@ -272,7 +275,7 @@
allocate(elem_proc_null(nproc_null))
! Filling up proc = 0 elements
nproc_null = 0
- do i = 1,nelmnts
+ do i = 0,nelmnts-1
if ( part(i) == 0 ) then
nproc_null = nproc_null + 1
elem_proc_null(nproc_null) = i
@@ -284,17 +287,20 @@
!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
-
+ !pgb: Solution , fault parellelization .
+ ipart=1
+ if (nproc > 1) then
+ do i = 1, nproc_null
+ part(elem_proc_null(i)) = ipart
+ if ( ipart == nproc-1 ) ipart = 0
+ ipart = ipart +1
+ end do
+ end if
+ deallocate(elem_proc_null)
endif
- call mesh2dual_ncommonnodes_fault(nelmnts, nnodes, nsize, sup_neighbour, &
- elmnts, xadj, adjncy, nnodes_elmnts, &
- nodes_elmnts, max_neighbour, 4, esize)
+! call mesh2dual_ncommonnodes_fault(nelmnts, nnodes, nsize, sup_neighbour, &
+! elmnts, xadj, adjncy, nnodes_elmnts, &
+! nodes_elmnts, max_neighbour, 4, esize)
! coupled elements
! ---------------
@@ -310,21 +316,58 @@
! 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
+!
+! =========== FAULT SHIELD double-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
+! ===================================================== !
+! Fault zone layer = the set of elements that contain at least one fault node
+ print *, "Fault zone layer :"
+
+! List of elements per node
+! nnodes_elmnts(i) = number of elements containing node #i (i=0:nnodes-1)
+! nodes_elmnts(nsize*i:nsize*i+nnodes_elmnts(i)-1) = index of elements (starting at 0) containing node #i
+! nsize = maximun number of elements in a node.
+! esize = nodes per element.
+
+ nnodes_elmnts(:) = 0
+ nodes_elmnts(:) = 0
+
+ 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
+
+ do iflt=1,size(faults)
+ do e=1,faults(iflt)%nspec
+ do k=1,4
+
+ inode = faults(iflt)%inodes1(k,e)-1 ! node index, starting at 0
+ k1 = nsize*inode
+ k2 = k1 + nnodes_elmnts(inode) -1
+ part( nodes_elmnts(k1:k2) ) = 0
+ inode = faults(iflt)%inodes2(k,e)-1
+ k1 = nsize*inode
+ k2 = k1 + nnodes_elmnts(inode) -1
+ part( nodes_elmnts(k1:k2) ) = 0
+
+ end do
+ end do
+ end do
+
nproc_null_final = count( part == 0 )
print *, nproc_null_final
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_movie_shakemap_AVS_DX_GMT.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_movie_shakemap_AVS_DX_GMT.f90 2011-06-28 12:00:53 UTC (rev 18663)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_movie_shakemap_AVS_DX_GMT.f90 2011-06-28 14:47:58 UTC (rev 18664)
@@ -568,11 +568,11 @@
max_field_current = + max_absol
! normalize field to [0:1]
- if( abs(max_field_current - min_field_current) > TINYVAL ) &
- field_display(:) = (field_display(:) - min_field_current) / (max_field_current - min_field_current)
+ ! if( abs(max_field_current - min_field_current) > TINYVAL ) &
+ ! field_display(:) = (field_display(:) - min_field_current) / (max_field_current - min_field_current)
! rescale to [-1,1]
- field_display(:) = 2.*field_display(:) - 1.
+ ! field_display(:) = 2.*field_display(:) - 1.
! apply threshold to normalized field
if(APPLY_THRESHOLD) &
@@ -588,10 +588,11 @@
endif
! map back to [0,1]
- field_display(:) = (field_display(:) + 1.) / 2.
+ ! field_display(:) = (field_display(:) + 1.) / 2.
+! Percy , Scale does not need to be rescaled , real units needs to be shown in the movie.
! map field to [0:255] for AVS color scale
- field_display(:) = 255. * field_display(:)
+! field_display(:) = 255. * field_display(:)
! apply scaling only if selected for shaking map
More information about the CIG-COMMITS
mailing list