[cig-commits] r18577 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db/src: . devel
percygalvez at geodynamics.org
percygalvez at geodynamics.org
Fri Jun 10 09:11:45 PDT 2011
Author: percygalvez
Date: 2011-06-10 09:11:44 -0700 (Fri, 10 Jun 2011)
New Revision: 18577
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/devel/fault_solver.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
splay_faults
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/devel/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/devel/fault_solver.f90 2011-06-10 16:07:25 UTC (rev 18576)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/devel/fault_solver.f90 2011-06-10 16:11:44 UTC (rev 18577)
@@ -89,15 +89,12 @@
!Number of time steps defined by the user : NTOUT
integer, save :: NTOUT,NSNAP
- integer, save :: SIMULATION_TYPE_DYN = 1
+ logical, save :: SIMULATION_TYPE_DYN = .false.
-
- integer , save :: size_Kelvin_Voigt
-
real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
- size_Kelvin_Voigt, SIMULATION_TYPE_DYN
+ SIMULATION_TYPE_DYN
contains
@@ -120,6 +117,8 @@
real(kind=CUSTOM_REAL) :: dt
integer :: iflt,ier,dummy_idfault
integer :: nbfaults
+ integer :: size_Kelvin_Voigt
+ integer :: SIMULATION_TYPE
character(len=256) :: filename
integer, parameter :: IIN_PAR =151
integer, parameter :: IIN_BIN =170
@@ -128,49 +127,61 @@
dummy_idfault = 0
- filename = prname(1:len_trim(prname))//'Kelvin_voigt_eta.bin'
- open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
- if( ier /= 0 ) stop 'Haven not found proc*_Kelvin_voigt_eta.bin'
- read(IIN_BIN) size_Kelvin_Voigt
- if (size_Kelvin_Voigt > 0) then
- allocate(Kelvin_Voigt_eta(size_Kelvin_Voigt))
- read(IIN_BIN) Kelvin_Voigt_eta
- endif
- Close(IIN_BIN)
-
open(unit=IIN_PAR,file='DATA/FAULT/Par_file_faults.in',status='old',iostat=ier)
if( ier /= 0 ) then
- write(6,*) 'Have not found Par_file_faults.in: assume no faults'
+ write(6,*) 'File Par_file_faults.in not found: assume no faults'
+ close(IIN_PAR)
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'
+ if( ier /= 0 ) then
+ write(6,*) 'File ',trim(filename),' not found. Abort'
+ stop
+ endif
! WARNING TO DO: should be an MPI abort
read(IIN_PAR,*) nbfaults
+ if (nbfaults==0) return
+! Reading etas of each fault
+ do iflt = 1,nbfaults
+ read(IIN_PAR,*) ! etas
+ enddo
+ read(IIN_PAR,*) SIMULATION_TYPE
+ if ( SIMULATION_TYPE /= 1 ) then
+ close(IIN_BIN)
+ close(IIN_PAR)
+ return
+ endif
+ SIMULATION_TYPE_DYN = .true.
+ read(IIN_PAR,*) NTOUT
+ read(IIN_PAR,*) NSNAP
+ read(IIN_PAR,*) V_HEALING
+ read(IIN_PAR,*) V_RUPT
+
+ read(IIN_BIN) nbfaults ! should be the same as in IIN_PAR
+ allocate( faults(nbfaults) )
do iflt=1,nbfaults
- read(IIN_PAR,*)
+ read(IIN_PAR,nml=BEGIN_FAULT,end=100)
+ call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
enddo
- read(IIN_PAR,*) SIMULATION_TYPE_DYN
- if ( SIMULATION_TYPE_DYN == 1 ) then
- read(IIN_PAR,*) NTOUT
- read(IIN_PAR,*) NSNAP
- read(IIN_PAR,*) V_HEALING
- read(IIN_PAR,*) V_RUPT
-
- read(IIN_BIN) nbfaults ! should be the same as in IIN_PAR
- allocate( faults(nbfaults) )
- 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)
- enddo
- endif
close(IIN_BIN)
close(IIN_PAR)
-
+
+ filename = prname(1:len_trim(prname))//'Kelvin_voigt_eta.bin'
+ open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ write(6,*) 'File ',trim(filename),' not found. Abort'
+ stop
+ endif
+ read(IIN_BIN) size_Kelvin_Voigt
+ if (size_Kelvin_Voigt > 0) then
+ allocate(Kelvin_Voigt_eta(size_Kelvin_Voigt))
+ read(IIN_BIN) Kelvin_Voigt_eta
+ endif
+ close(IIN_BIN)
return
100 stop 'Did not find BEGIN_FAULT block #'
! WARNING TO DO: should be an MPI abort
@@ -203,7 +214,6 @@
read(IIN_BIN) bc%nspec,bc%nglob
if (bc%nspec==0) return
-
allocate( bc%ibulk1(bc%nglob) )
allocate( bc%ibulk2(bc%nglob) )
allocate( ibool1(NGLLSQUARE,bc%nspec) )
@@ -450,15 +460,14 @@
real(kind=CUSTOM_REAL), intent(in) :: V(:,:),D(:,:)
integer,intent(in) :: iflt
-
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength
real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: t1,t2,tnorm,tnew
real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: dD,dV,dA
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, Vnorm, Vnorm_old
real(kind=CUSTOM_REAL) :: half_dt
-! integer :: k
-
+! integer :: k
+
half_dt = 0.5e0_CUSTOM_REAL*bc%dt
theta_old = bc%swf%theta
Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
@@ -508,8 +517,10 @@
! Solve for shear stress
tnorm = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
+ tnorm = max(tnorm,1e0_CUSTOM_REAL)
t1 = T(1,:)/tnorm
t2 = T(2,:)/tnorm
+
tnew = min(tnorm,strength)
T(1,:) = tnew * t1
T(2,:) = tnew * t2
@@ -910,7 +921,6 @@
if (vold(i)<=V_RUPT .and. vnew(i)>=V_RUPT) dataXZ%tRUP(i)= time-dt*(vnew(i)-V_RUPT)/(vnew(i)-vold(i))
endif
enddo
-
! To do : add stress criteria (firs time strength is reached).
@@ -930,34 +940,28 @@
write(filename,"('OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
- open(unit=IOUT, file= trim(filename), status='replace', form='formatted',action='write')
-! open(unit=IOUT, file= trim(filename), status='replace', form='unformatted')
+! open(unit=IOUT, file= trim(filename), status='replace', form='formatted',action='write')
! NOTE : It had to be adopted formatted output to avoid conflicts readings with different
! compilers.
- write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
+! write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
+
+ open(unit=IOUT, file= trim(filename), status='replace', form='unformatted',action='write')
-
-! WARNING: for the case of multiple faults the filename must contain a fault identifier
-! (a separate snapshot file for each fault)
-! write(filename,"('OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
-!
-! open(unit=IOUT, file= trim(filename), status='replace', form='unformatted')
-
-! write(IOUT) dataXZ%xcoord
-! write(IOUT) dataXZ%ycoord
-! write(IOUT) dataXZ%zcoord
-! write(IOUT) dataXZ%d1
-! write(IOUT) dataXZ%d2
-! write(IOUT) dataXZ%v1
-! write(IOUT) dataXZ%v2
-! write(IOUT) dataXZ%t1
-! write(IOUT) dataXZ%t2
-! write(IOUT) dataXZ%t3
-! write(IOUT) dataXZ%sta
-! write(IOUT) dataXZ%stg
-! write(IOUT) dataXZ%tRUP
-! write(IOUT) dataXZ%tPZ
+ write(IOUT) dataXZ%xcoord
+ write(IOUT) dataXZ%ycoord
+ write(IOUT) dataXZ%zcoord
+ write(IOUT) dataXZ%d1
+ write(IOUT) dataXZ%d2
+ write(IOUT) dataXZ%v1
+ write(IOUT) dataXZ%v2
+ write(IOUT) dataXZ%t1
+ write(IOUT) dataXZ%t2
+ write(IOUT) dataXZ%t3
+ write(IOUT) dataXZ%sta
+ write(IOUT) dataXZ%stg
+ write(IOUT) dataXZ%tRUP
+ write(IOUT) dataXZ%tPZ
close(IOUT)
end subroutine write_dataXZ
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 2011-06-10 16:07:25 UTC (rev 18576)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2011-06-10 16:11:44 UTC (rev 18577)
@@ -86,28 +86,26 @@
integer :: nb,i,iflt,ier,nspec,dummy_node
integer, parameter :: IIN = 100
- real(kind=CUSTOM_REAL) :: eta
! read fault input file
nb = 0
open(unit=IIN,file='DATA/FAULT/Par_file_faults.in',status='old',action='read',iostat=ier)
if (ier==0) then
read(IIN,*) nb
- read(IIN,*) eta
else
- write(6,*) 'File Par_file_faults.in does not exist '
- return
+ write(6,*) 'No faults in the domain'
+ write(6,*) 'Par_file_faults.in does not exist '
+ close(IIN)
end if
- close(IIN)
ANY_FAULT = (nb>0)
+ if (.not. ANY_FAULT) return
- if (.not. ANY_FAULT) return
-
allocate(fault_db(nb))
do i=1,nb
- fault_db(i)%eta = eta
+ read(IIN,*) fault_db(i)%eta
enddo
+ close(IIN)
! read fault database file
open(unit=IIN,file=prname(1:len_trim(prname))//'Database_fault', &
@@ -178,7 +176,7 @@
do iflt=1,size(fault_db)
! TO DO : Percy, iface is computed with fault open , we have to change this and close *_dummy
-! at this point.
+! at this point.
call setup_iface(fault_db(iflt),nnodes_ext_mesh,nodes_coords_ext_mesh,nspec,nglob,ibool)
! saving gll indices for each fault face, needed for ibulks
@@ -616,16 +614,16 @@
if (.not.ANY_FAULT) return
! opening Kelvin_voig_eta.bin (Necessary for all processors , if number of fault elements = 0
! then the file will be empty with
- filename = prname(1:len_trim(prname))//'Kelvin_voigt_eta.bin'
+ filename = prname(1:len_trim(prname))//'Kelvin_voigt_eta.bin'
open(unit=IOUT,file=trim(filename),status='unknown',action='write',form='unformatted',iostat=ier)
if( ier /= 0 ) then
- write(IOUT,*) 'error opening file ',trim(filename)
+ write(IOUT) 'error opening file ',trim(filename)
stop
endif
! saves mesh file proc***_Kelvin_voigt_eta.bin
if (allocated(Kelvin_Voigt_eta)) then
- write(IOUT) size(Kelvin_Voigt_eta)
+ size_Kelvin_Voigt = size(Kelvin_Voigt_eta)
else
size_Kelvin_Voigt = 0
endif
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 2011-06-10 16:07:25 UTC (rev 18576)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2011-06-10 16:11:44 UTC (rev 18577)
@@ -129,7 +129,8 @@
open(unit=IIN_PAR,file='DATA/FAULT/Par_file_faults.in',status='old',iostat=ier)
if( ier /= 0 ) then
- write(6,*) 'File Par_file_faults.in not found: assume no faults'
+ write(6,*) 'File Par_file_faults.in not found: assume no faults'
+ close(IIN_PAR)
return
endif
@@ -144,7 +145,10 @@
read(IIN_PAR,*) nbfaults
if (nbfaults==0) return
- read(IIN_PAR,*) ! eta
+! Reading etas of each fault
+ do iflt = 1,nbfaults
+ read(IIN_PAR,*) ! etas
+ enddo
read(IIN_PAR,*) SIMULATION_TYPE
if ( SIMULATION_TYPE /= 1 ) then
close(IIN_BIN)
@@ -165,7 +169,7 @@
enddo
close(IIN_BIN)
close(IIN_PAR)
-
+
filename = prname(1:len_trim(prname))//'Kelvin_voigt_eta.bin'
open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) then
@@ -178,7 +182,6 @@
read(IIN_BIN) Kelvin_Voigt_eta
endif
close(IIN_BIN)
-
return
100 stop 'Did not find BEGIN_FAULT block #'
! WARNING TO DO: should be an MPI abort
@@ -211,7 +214,6 @@
read(IIN_BIN) bc%nspec,bc%nglob
if (bc%nspec==0) return
-
allocate( bc%ibulk1(bc%nglob) )
allocate( bc%ibulk2(bc%nglob) )
allocate( ibool1(NGLLSQUARE,bc%nspec) )
@@ -458,15 +460,14 @@
real(kind=CUSTOM_REAL), intent(in) :: V(:,:),D(:,:)
integer,intent(in) :: iflt
-
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength
real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: t1,t2,tnorm,tnew
real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: dD,dV,dA
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, Vnorm, Vnorm_old
real(kind=CUSTOM_REAL) :: half_dt
-! integer :: k
-
+! integer :: k
+
half_dt = 0.5e0_CUSTOM_REAL*bc%dt
theta_old = bc%swf%theta
Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
@@ -516,8 +517,10 @@
! Solve for shear stress
tnorm = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
+ tnorm = max(tnorm,1e0_CUSTOM_REAL)
t1 = T(1,:)/tnorm
t2 = T(2,:)/tnorm
+
tnew = min(tnorm,strength)
T(1,:) = tnew * t1
T(2,:) = tnew * t2
@@ -918,7 +921,6 @@
if (vold(i)<=V_RUPT .and. vnew(i)>=V_RUPT) dataXZ%tRUP(i)= time-dt*(vnew(i)-V_RUPT)/(vnew(i)-vold(i))
endif
enddo
-
! To do : add stress criteria (firs time strength is reached).
@@ -938,34 +940,28 @@
write(filename,"('OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
- open(unit=IOUT, file= trim(filename), status='replace', form='formatted',action='write')
-! open(unit=IOUT, file= trim(filename), status='replace', form='unformatted')
+! open(unit=IOUT, file= trim(filename), status='replace', form='formatted',action='write')
! NOTE : It had to be adopted formatted output to avoid conflicts readings with different
! compilers.
- write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
-
-
-! WARNING: for the case of multiple faults the filename must contain a fault identifier
-! (a separate snapshot file for each fault)
-! write(filename,"('OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
-!
-! open(unit=IOUT, file= trim(filename), status='replace', form='unformatted')
+! write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
-! write(IOUT) dataXZ%xcoord
-! write(IOUT) dataXZ%ycoord
-! write(IOUT) dataXZ%zcoord
-! write(IOUT) dataXZ%d1
-! write(IOUT) dataXZ%d2
-! write(IOUT) dataXZ%v1
-! write(IOUT) dataXZ%v2
-! write(IOUT) dataXZ%t1
-! write(IOUT) dataXZ%t2
-! write(IOUT) dataXZ%t3
-! write(IOUT) dataXZ%sta
-! write(IOUT) dataXZ%stg
-! write(IOUT) dataXZ%tRUP
-! write(IOUT) dataXZ%tPZ
+ open(unit=IOUT, file= trim(filename), status='replace', form='unformatted',action='write')
+
+ write(IOUT) dataXZ%xcoord
+ write(IOUT) dataXZ%ycoord
+ write(IOUT) dataXZ%zcoord
+ write(IOUT) dataXZ%d1
+ write(IOUT) dataXZ%d2
+ write(IOUT) dataXZ%v1
+ write(IOUT) dataXZ%v2
+ write(IOUT) dataXZ%t1
+ write(IOUT) dataXZ%t2
+ write(IOUT) dataXZ%t3
+ write(IOUT) dataXZ%sta
+ write(IOUT) dataXZ%stg
+ write(IOUT) dataXZ%tRUP
+ write(IOUT) dataXZ%tPZ
close(IOUT)
end subroutine write_dataXZ
More information about the CIG-COMMITS
mailing list