[cig-commits] r18349 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: decompose_mesh_SCOTCH src
ampuero at geodynamics.org
ampuero at geodynamics.org
Wed May 11 18:27:16 PDT 2011
Author: ampuero
Date: 2011-05-11 18:27:16 -0700 (Wed, 11 May 2011)
New Revision: 18349
Modified:
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_kinematic.f90
Log:
tried fixing write_fault_database changing glob2loc* to pointers
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-12 00:23:10 UTC (rev 18348)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2011-05-12 01:27:16 UTC (rev 18349)
@@ -12,10 +12,7 @@
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_DOUBLE
- double precision, parameter :: FAULT_GAP_TOLERANCE = 1e0_CUSTOM_REAL
+ double precision, parameter :: FAULT_GAP_TOLERANCE = 1e0 ! JPA: are you sure 1 meter is small enough?
public :: read_fault_files, fault_collect_elements, close_faults, write_fault_database
@@ -58,53 +55,47 @@
character(len=5) :: NTchar
integer :: e,ier,nspec_side1,nspec_side2
- write(NTchar,'(I5)') ifault
- NTchar = adjustl(NTchar)
+ 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)
+ 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.'
- 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)
- !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))
- !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
- else
+ 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
@@ -170,8 +161,8 @@
dist = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
dist = sqrt(dist)
- if ((0.0_CUSTOM_REAL < dist) .and. (dist <= FAULT_GAP_TOLERANCE)) then
- xyz = (xyz_1 + xyz_2)*0.5_CUSTOM_REAL
+ if (dist <= FAULT_GAP_TOLERANCE) then
+ xyz = (xyz_1 + xyz_2)*0.5d0
nodes_coords(:,iglob2) = xyz
nodes_coords(:,iglob1) = xyz
found_it = .true.
@@ -332,9 +323,9 @@
integer(long), intent(in) :: nelmnts
integer, dimension(1:nelmnts), intent(in) :: part
integer, dimension(0:nelmnts-1), intent(in) :: glob2loc_elmnts
- integer, dimension(:), intent(in) :: glob2loc_nodes_nparts
- integer, dimension(:), intent(in) :: glob2loc_nodes_parts
- integer, dimension(:), intent(in) :: glob2loc_nodes
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
integer :: i,j,k,iflt,e
integer :: nspec_fault_1,nspec_fault_2
@@ -352,14 +343,15 @@
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) == iproc) then
+ if (part(e) == 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
@@ -372,6 +364,7 @@
end if
enddo
+ ! export fault element data, side 2
do i=1,faults(iflt)%nspec
e = faults(iflt)%ispec2(i)
if(part(e) == iproc) then
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 2011-05-12 00:23:10 UTC (rev 18348)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2011-05-12 01:27:16 UTC (rev 18349)
@@ -57,7 +57,7 @@
integer :: nspec,nglob
real(kind=CUSTOM_REAL) :: dt
real(kind=CUSTOM_REAL), dimension(:), pointer :: B,invM1,invM2,Z
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T,slip,slip_rate,coord
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T,D,V,coord
real(kind=CUSTOM_REAL), dimension(:,:,:), pointer :: R
integer, dimension(:), pointer :: ibulk1, ibulk2
type(dataT_type) :: dataT
@@ -224,12 +224,12 @@
bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
allocate(bc%T(3,bc%nglob))
- allocate(bc%slip(3,bc%nglob))
- allocate(bc%slip_rate(3,bc%nglob))
+ allocate(bc%D(3,bc%nglob))
+ allocate(bc%V(3,bc%nglob))
bc%T = 0e0_CUSTOM_REAL
- bc%slip = 0e0_CUSTOM_REAL
- bc%slip_rate = 0e0_CUSTOM_REAL
-! Dt between two loaded slip_rates
+ bc%D = 0e0_CUSTOM_REAL
+ bc%V = 0e0_CUSTOM_REAL
+! Dt between two loaded slip rates
read(IIN_PAR,nml=KINPAR)
bc%kin_dt = kindt
@@ -369,28 +369,28 @@
endif
! Kinematic velocity_rate
-! bc%slip_rate : Imposed apriori and read from slip rate snapshots (from time reversal)
+! bc%V : Imposed apriori and read from slip rate snapshots (from time reversal)
! linear interpolate between consecutive kinematic time steps.
-! slip_rate will be given each time step.
- bc%slip_rate(1,:) = ( (t2 - time)*bc%v_kin_t1(1,:) + (time - t1)*bc%v_kin_t2(1,:) )/ bc%kin_dt
- bc%slip_rate(2,:) = ( (t2 - time)*bc%v_kin_t1(2,:) + (time - t1)*bc%v_kin_t2(2,:) )/ bc%kin_dt
+! V will be given each time step.
+ bc%V(1,:) = ( (t2 - time)*bc%v_kin_t1(1,:) + (time - t1)*bc%v_kin_t2(1,:) )/ bc%kin_dt
+ bc%V(2,:) = ( (t2 - time)*bc%v_kin_t1(2,:) + (time - t1)*bc%v_kin_t2(2,:) )/ bc%kin_dt
!dV_free = dV_predictor + (dt/2)*dA_free
dV_free(1,:) = dV(1,:)+half_dt*dA(1,:)
dV_free(2,:) = dV(2,:)+half_dt*dA(2,:)
dV_free(3,:) = dV(3,:)+half_dt*dA(3,:)
-! T = Z*( dV_free - V_slip_rate) , V_slip_rate known apriori as input.
+! T = Z*( dV_free - V) , V known apriori as input.
! CONVENTION : T(ibulk1)=T=-T(ibulk2)
- T(1,:) = bc%Z * ( dV_free(1,:) -bc%slip_rate(1,:) )
- T(2,:) = bc%Z * ( dV_free(2,:) -bc%slip_rate(2,:) )
+ T(1,:) = bc%Z * ( dV_free(1,:) -bc%V(1,:) )
+ T(2,:) = bc%Z * ( dV_free(2,:) -bc%V(2,:) )
T(3,:) = bc%Z * ( dV_free(3,:) )
! Save tractions
bc%T = T
! Update slip in fault frame
- bc%slip = dD
+ bc%D = dD
! Rotate tractions back to (x,y,z) frame
T = rotate(bc,T,-1)
@@ -405,7 +405,7 @@
MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
!-- intermediate storage of outputs --
- call store_dataT(bc%dataT,bc%slip,bc%slip_rate,bc%T,it)
+ call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
!-- OUTPUTS --
! write dataT every NTOUT time steps or at the end of simulation
@@ -648,10 +648,10 @@
! coord : Receivers coordinates
! npoin : number of Receivers.
! nglob : number of gll points along the fault.
-! dataXZ : Velocity slip_rate .
+! dataXZ : slip rate .
! iflt : number of faults.
-! OUTPUT v : slip_rate on receivers.
+! OUTPUT v : slip rate on receivers.
subroutine load_vslip_snapshots(dataXZ,itime,nglob,iflt)
More information about the CIG-COMMITS
mailing list