[cig-commits] r20747 - in seismo/3D/FAULT_SOURCE: . branches/new_fault_db/CUBIT branches/new_fault_db/decompose_mesh_SCOTCH branches/new_fault_db/src
ampuero at geodynamics.org
ampuero at geodynamics.org
Wed Sep 19 17:01:30 PDT 2012
Author: ampuero
Date: 2012-09-19 17:01:29 -0700 (Wed, 19 Sep 2012)
New Revision: 20747
Modified:
seismo/3D/FAULT_SOURCE/TO_DO
seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/functions.py
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_solver.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90
Log:
fixed fault_repartition_parallel for triple junctions; some clean up
Modified: seismo/3D/FAULT_SOURCE/TO_DO
===================================================================
--- seismo/3D/FAULT_SOURCE/TO_DO 2012-09-19 23:16:22 UTC (rev 20746)
+++ seismo/3D/FAULT_SOURCE/TO_DO 2012-09-20 00:01:29 UTC (rev 20747)
@@ -1,22 +1,15 @@
Development of dynamic and kinematic earthquake rupture in SPECFEM3D-SESAME.
Prioritized to-do list:
-+ implement rate-and-state friction, with option for several versions of the friction law
++ parallelized fault:
+ - debug the issue with nodes shared by three processors
+ - in fault_solver.f90: parallelize dataT outputs
++ rate-and-state friction:
+ - debug
+ - option for several versions of the friction law
+ verify consistency of fault edge nodes in kinematic solver:
currently, these nodes are not split but they are included in the fault database
-+ parallelize the fault:
- - in fault_scotch.f90:
- . reorder the fault element lists so that ispec1 and ispec2
- are facing each other: for each face ispec1 search the
- face ispec2 with nearest center point
- . for each fault, for each pair (ispec1,ispec2) put both
- elements in the processor with lowest rank among the two
- . repeat the loop above to handle triple junctions
- - in fault_solver.f90:
- . assemble B and the fault normal vectors during initialization
- . make sure internal forces are assembled across processors before fault solver
- No need to re-assemble after fault solver because B is already assembled
- . parallelize outputs
+
+ simplify mesh generation with faults:
- eliminate fault edge nodes?
Currently fault edge nodes must be defined as closed and non-split in CUBIT,
@@ -24,6 +17,9 @@
- split nodes in Scotch/SESAME instead of in CUBIT ?
- develop a recipe for mesh coarsening in CUBIT
++ cubit interface:
+ - consolidate save_fault_nodes_elements.py and functions.py
+
+ add snapshot outputs to kinematic solver (e.g. to export fault stresses)
+ factor common subroutines of fault solvers into a lower level module
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/functions.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/functions.py 2012-09-19 23:16:22 UTC (rev 20746)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/functions.py 2012-09-20 00:01:29 UTC (rev 20747)
@@ -101,14 +101,10 @@
fsideu = 'MESH_tohoku_break_double_surface/fault_sideu.dat'
fsided = 'MESH_tohoku_break_double_surface/fault_sided.dat'
+delta = 0.5 # WARNING: Make sure that this delta is smaller than FAULT_GAP_TOLERANCE
+ # defined in decompose_mesh_SCOTH/fault_scotch.f90 (usually FAULT_GAP_TOLERANCE=1.0d0)
+ # and larger than SMALLVAL_TOL defined in constants.h (usually SMALLVAL_TOL=1.d-10*dabs(UTM_X_MAX - UTM_X_MIN))
-delta = 0.5 # Make sure that this delta is in coorcondance with
- # FAULT_GAP_TOLERANCE which is normally set up to 1.0d0 but can be changed
- # in decompose_mesh_SCOTH/fault_scotch.f90, and higher
- # than constants.h/SMALLVAL_TOL used by the routine get_global.f90.
- # SMALLVALTOL = 1.d-10 * dabs(UTM_X_MAX - UTM_X_MIN)
-
-
def m2km():
name_in = 'MESH/nodes_coords_file'
name_out = 'nodes_coords_file_m'
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-19 23:16:22 UTC (rev 20746)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-09-20 00:01:29 UTC (rev 20747)
@@ -395,9 +395,7 @@
call save_nodes_coords(nodes_coords,nnodes)
call close_faults(nodes_coords,elmnts,nspec,nnodes,esize)
!SNS
- if(PARALLEL_FAULT) then
- call reorder_fault_elements(nodes_coords,elmnts,nnodes,esize,nspec)
- endif
+ if (PARALLEL_FAULT) call reorder_fault_elements(nodes_coords,nnodes)
end if
end subroutine read_mesh_files
@@ -538,7 +536,6 @@
if(.not. PARALLEL_FAULT) then
! 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, nnodes, elmnts, nsize, nparts, part, esize)
endif
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-19 23:16:22 UTC (rev 20746)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2012-09-20 00:01:29 UTC (rev 20747)
@@ -203,141 +203,82 @@
end subroutine close_fault_single
!===================================================================================================
- subroutine reorder_fault_elements(nodes_coords,elmnts,nnodes,esize,nelmnts)
- integer ,intent(in) :: nnodes, esize
- integer(long), intent(in) :: nelmnts
- integer, dimension(esize,nelmnts), intent(in) :: elmnts
- double precision,dimension(3,nnodes), intent(inout) :: nodes_coords
+ subroutine reorder_fault_elements(nodes_coords,nnodes)
- integer :: iflt
+ integer, intent(in) :: nnodes
+ double precision,dimension(3,nnodes), intent(in) :: nodes_coords
+
+ integer :: i
- do iflt=1,size(faults)
- call reorder_fault_elements_single(faults(iflt),&
- elmnts,nodes_coords,nnodes,esize,nelmnts,faults(iflt)%nspec)
+ do i=1,size(faults)
+ call reorder_fault_elements_single(faults(i),nodes_coords,nnodes)
enddo
end subroutine reorder_fault_elements
! ---------------------------------------------------------------------------------------------------
- subroutine reorder_fault_elements_single(f,elmnts,nodes_coords,nnodes,esize,nelmnts,nspec)
+ subroutine reorder_fault_elements_single(f,nodes_coords,nnodes)
- integer ,intent(in) :: nnodes, esize
- integer(long), intent(in) :: nelmnts
- integer, dimension(esize,nelmnts), intent(in) :: elmnts
- double precision,dimension(3,nnodes), intent(inout) :: nodes_coords
type(fault_type), intent(inout) :: f
+ integer, intent(in) :: nnodes
+ double precision, dimension(3,nnodes), intent(in) :: nodes_coords
- double precision, dimension(3,4) :: xyz
- double precision, dimension(:,:), allocatable :: xyz_c
+ double precision, dimension(3,f%nspec) :: xyz_c
+ integer :: k,e,iglob
+ integer, dimension(f%nspec) :: new_index_list
- double precision :: dist
- integer :: iglob1, iglob2, i, j, k
- logical :: found_it
-
- integer :: e,iglob
-
- integer :: nspec
-
- integer, dimension(nspec) :: new_index_list
-
- allocate(xyz_c(3,f%nspec))
- xyz_c=0._CUSTOM_REAL
-
! compute element-face centers for fault side 1
+ xyz_c = 0d0
do e=1,f%nspec
do k=1,4
iglob = f%inodes1(k,e)
- xyz(:,k) = nodes_coords(:,iglob)
+ xyz_c(:,e) = xyz_c(:,e) + nodes_coords(:,iglob)
enddo
- xyz_c(1,e) = sum(xyz(1,:))/4
- xyz_c(2,e) = sum(xyz(2,:))/4
- xyz_c(3,e) = sum(xyz(3,:))/4
enddo
+ xyz_c = xyz_c / 4d0
! reorder
call lex_order(xyz_c,new_index_list,nnodes,f%nspec)
f%ispec1 = f%ispec1(new_index_list)
f%inodes1 = f%inodes1(:,new_index_list)
! repeat for fault side 2
+ xyz_c = 0d0
do e=1,f%nspec
do k=1,4
iglob = f%inodes2(k,e)
- xyz(:,k) = nodes_coords(:,iglob)
+ xyz_c(:,e) = xyz_c(:,e) + nodes_coords(:,iglob)
enddo
- xyz_c(1,e) = sum(xyz(1,:))/4
- xyz_c(2,e) = sum(xyz(2,:))/4
- xyz_c(3,e) = sum(xyz(3,:))/4
enddo
- ! reorder
+ xyz_c = xyz_c / 4d0
call lex_order(xyz_c,new_index_list,nnodes,f%nspec)
f%ispec2 = f%ispec2(new_index_list)
f%inodes2 = f%inodes2(:,new_index_list)
-
- do e=1,f%nspec
- do k=1,4
- iglob = f%inodes1(k,e)
- xyz(:,k) = nodes_coords(:,iglob)
- enddo
- xyz_c(1,e) = sum(xyz(1,:))/4
- xyz_c(2,e) = sum(xyz(2,:))/4
- xyz_c(3,e) = sum(xyz(3,:))/4
- enddo
-
-
- do e=1,f%nspec
- do k=1,4
- iglob = f%inodes2(k,e)
- xyz(:,k) = nodes_coords(:,iglob)
- enddo
- xyz_c(1,e) = sum(xyz(1,:))/4
- xyz_c(2,e) = sum(xyz(2,:))/4
- xyz_c(3,e) = sum(xyz(3,:))/4
- enddo
-
-
- deallocate(xyz_c)
-
end subroutine reorder_fault_elements_single
! ---------------------------------------------------------------------------------------------------
subroutine lex_order(xyz_c,loc,nnodes,nspec)
- integer ,intent(in) :: nnodes,nspec
- double precision, dimension(3,nspec) :: xyz_c
+ integer, intent(in) :: nnodes,nspec,loc(nspec)
+ double precision, intent(in) :: xyz_c(3,nspec)
- integer nglob !,nspec
- integer iglob(nnodes),loc(nspec)
- logical ifseg(nspec)
- double precision xp(nspec),yp(nspec),zp(nspec)
- double precision UTM_X_MIN,UTM_X_MAX
-
- integer ispec,i,j,ier
- integer ieoff,ilocnum,nseg,ioff,iseg,ig
-
- integer, dimension(:), allocatable :: ind,ninseg,iwork
- double precision, dimension(:), allocatable :: work
-
+ double precision, dimension(nspec) :: work,xp,yp,zp
+ integer, dimension(nspec) :: ind,ninseg,iwork
+ logical :: ifseg(nspec)
+ double precision :: UTM_X_MIN,UTM_X_MAX
+ integer :: ispec,i,j
+ integer :: ieoff,ilocnum,nseg,ioff,iseg
double precision :: SMALLVALTOL
-
xp=xyz_c(1,:)
yp=xyz_c(2,:)
zp=xyz_c(3,:)
+ ! define geometrical tolerance based upon typical size of the model
UTM_X_MAX = maxval(xp)
UTM_X_MIN = minval(xp)
- ! define geometrical tolerance based upon typical size of the model
SMALLVALTOL = 1.d-10 * dabs(UTM_X_MAX - UTM_X_MIN)
- ! dynamically allocate arrays
- allocate(ind(nspec), &
- ninseg(nspec), &
- iwork(nspec), &
- work(nspec),stat=ier)
- if( ier /= 0 ) stop 'error allocating arrays'
-
-
! establish initial pointers
do ispec=1,nspec
ieoff=1*(ispec-1)
@@ -396,12 +337,6 @@
enddo
enddo
- ! deallocate arrays
- deallocate(ind)
- deallocate(ninseg)
- deallocate(iwork)
- deallocate(work)
-
end subroutine lex_order
!===================================================================================================
@@ -409,11 +344,12 @@
! Repartitioning : two coupled faultside1/side2 elements are transfered to the same partition
!--------------------------------------------------
- Subroutine fault_repartition (nelmnts, nnodes, elmnts, nsize, &
+ subroutine fault_repartition (nelmnts, nnodes, elmnts, nsize, &
nproc, part, esize)
-! 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.
+! part(e) = index of the processor to which element #e is assigned.
+! Fault elements and neighbors are assigned to the same processor.
+! Part, once modified, will be input for write_partition_database.
!INPUTS
integer(long), intent(in) :: nelmnts,nsize
@@ -424,7 +360,6 @@
integer, dimension(0:nnodes-1) :: nnodes_elmnts
integer, dimension(0:nsize*nnodes-1) :: nodes_elmnts
- integer :: max_neighbour
integer :: i,j, ipart,nproc_null,nproc_null_final
integer :: k1, k2, k,e,iflt,inode
@@ -499,46 +434,29 @@
end subroutine fault_repartition
! ---------------------------------------------------------------------------------------------------
- subroutine fault_repartition_parallel (nelmnts, nnodes, elmnts, nsize, &
- nproc, part, esize)
+ subroutine fault_repartition_parallel (nelmnts, part)
+ integer(long), intent(in) :: nelmnts
+ integer, dimension(0:nelmnts-1), intent(inout) :: part
+ integer :: i,iflt,e,e1,e2,proc1,proc2
-!INPUTS
- integer(long), intent(in) :: nelmnts,nsize
- integer, intent(in) :: nnodes, nproc, esize
- integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
-!OUTPUTS :
- integer, dimension(0:nelmnts-1), intent(inout) :: part
+!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
- integer, dimension(0:nnodes-1) :: nnodes_elmnts
- integer, dimension(0:nsize*nnodes-1) :: nodes_elmnts
- integer :: max_neighbour
-
- integer :: i,j, ipart,nproc_null,nproc_null_final
- integer :: k1, k2, k,e,iflt,inode
- integer, dimension(:), allocatable :: elem_proc_null
-
-
- integer :: e1,e2,proc1,proc2
-
- do iflt=1,size(faults)
- do e=1,faults(iflt)%nspec
- e1 = faults(iflt)%ispec1(e) - 1
- e2 = faults(iflt)%ispec2(e) - 1
- proc1 = part(e1)
- proc2 = part(e2)
- part(e1) = min(proc1,proc2)
- part(e2) = part(e1)
+ do i=1,2 !JPA do it twice, to handle triple junctions (intersections between two faults)
+ do iflt=1,size(faults)
+ do e=1,faults(iflt)%nspec
+ e1 = faults(iflt)%ispec1(e) - 1
+ e2 = faults(iflt)%ispec2(e) - 1
+ proc1 = part(e1)
+ proc2 = part(e2)
+ part(e1) = min(proc1,proc2)
+ part(e2) = part(e1)
+ end do
end do
- end do
-
- do iflt=1,size(faults)
- do e=1,faults(iflt)%nspec
- e1 = faults(iflt)%ispec1(e) - 1
- e2 = faults(iflt)%ispec2(e) - 1
- enddo
- enddo
+ enddo
end subroutine fault_repartition_parallel
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-09-19 23:16:22 UTC (rev 20746)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-09-20 00:01:29 UTC (rev 20747)
@@ -164,7 +164,6 @@
return
endif
- dt = real(DTglobal)
filename = prname(1:len_trim(prname))//'fault_db.bin'
open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) then
@@ -193,6 +192,7 @@
read(IIN_BIN) nbfaults ! should be the same as in IIN_PAR
allocate( faults(nbfaults) )
+ dt = real(DTglobal)
do iflt=1,nbfaults
read(IIN_PAR,nml=BEGIN_FAULT,end=100)
call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt,accel)
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2012-09-19 23:16:22 UTC (rev 20746)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2012-09-20 00:01:29 UTC (rev 20747)
@@ -97,7 +97,7 @@
real(kind=CUSTOM_REAL) :: dt
integer :: iflt,ier,dummy_idfault
integer :: nbfaults
- integer :: SIMULATION_TYPE
+ integer :: SIMULATION_TYPE
character(len=256) :: filename
integer, parameter :: IIN_PAR =151
integer, parameter :: IIN_BIN =170
@@ -113,7 +113,6 @@
return
endif
- dt = real(DTglobal)
filename = prname(1:len_trim(prname))//'fault_db.bin'
open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) stop 'Have not found proc*_fault_db.bin'
@@ -131,6 +130,7 @@
read(IIN_PAR,*) DUMMY
read(IIN_BIN) nbfaults ! should be the same as in IIN_PAR
allocate( faults(nbfaults) )
+ dt = real(DTglobal)
do iflt=1,nbfaults
read(IIN_PAR,nml=BEGIN_FAULT,end=100)
call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
More information about the CIG-COMMITS
mailing list