[cig-commits] r19682 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: . src
ampuero at geodynamics.org
ampuero at geodynamics.org
Sat Feb 25 08:39:12 PST 2012
Author: ampuero
Date: 2012-02-25 08:39:12 -0800 (Sat, 25 Feb 2012)
New Revision: 19682
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT
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:
cosmetic clean up
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT 2012-02-25 01:39:11 UTC (rev 19681)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT 2012-02-25 16:39:12 UTC (rev 19682)
@@ -135,6 +135,8 @@
nummaterial_velocity_file
The CUBIT graphics window should show the mesh (e.g. EXAMPLES/splay_faults/splay_faults.jpg)
+ You can zoom, pan or rotate using mouse gestures. To look at the block of your interest,
+ set "visibility off" in the other blocks.
To run CUBIT from different working directories it is convenient to
include the path to SPECFEM3D/CUBIT in an environment variable PYTHONPATH
@@ -426,12 +428,14 @@
VI. RUNNING ON FRAM (THE CALTECH GPS CLUSTER)
---------------------------------------------
+Getting started with Fram:
+ http://hpc.caltech.edu/hel/citerrafram-new-user-guide/
+ http://citerra.gps.caltech.edu/Fram_User_Meeting.pdf
+
Before step I.3 add this to you ~/.bash_profile file:
-
# Load modules:
module load intel/intel-12
module load intel/impi
-
# Set huge stack size:
ulimit -s unlimited
@@ -445,5 +449,8 @@
SCOTCH_LIBS = -L/home/surendra/decompose_mesh_SCOTCH/scotch_5.1/lib/ -lscotch -lscotcherr
Submitting a job:
- qsub -l nodes=48 -l walltime=30:00 -N simulation_name -m bae -M your at email.address go_mesher
+ qsub -l nodes=48 -l walltime=30:00 go_mesher
+Other useful options to qsub are:
+ -N simulation_name
+ -m bae -M your at email.address
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-02-25 01:39:11 UTC (rev 19681)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-02-25 16:39:12 UTC (rev 19682)
@@ -109,7 +109,7 @@
! Minv inverse mass matrix
! dt global time step
!
- subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt)
+subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt)
character(len=256), intent(in) :: prname ! 'proc***'
real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
@@ -197,7 +197,7 @@
!---------------------------------------------------------------------
- subroutine init_one_fault_TPV16(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
+subroutine init_one_fault_TPV16(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
type(bc_dynflt_type), intent(inout) :: bc
real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
@@ -416,7 +416,7 @@
!---------------------------------------------------------------------
- subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
+subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
type(bc_dynflt_type), intent(inout) :: bc
real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
@@ -554,10 +554,10 @@
call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
call init_dataXZ(bc%dataXZ,bc,bc%nglob)
- end subroutine init_one_fault
+end subroutine init_one_fault
!---------------------------------------------------------------------
- subroutine compute_R(R,nglob,nx,ny,nz)
+subroutine compute_R(R,nglob,nx,ny,nz)
integer :: nglob
real(kind=CUSTOM_REAL), intent(out) :: R(3,3,nglob)
@@ -590,11 +590,11 @@
R(3,2,:)=ny
R(3,3,:)=nz
- end subroutine compute_R
+end subroutine compute_R
!---------------------------------------------------------------------
! adds a value to a fault parameter inside an area with prescribed shape
- subroutine init_2d_distribution(a,coord,iin,n)
+subroutine init_2d_distribution(a,coord,iin,n)
real(kind=CUSTOM_REAL), intent(inout) :: a(:)
real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
@@ -649,10 +649,10 @@
where (b /= 0e0_CUSTOM_REAL) a = b
enddo
- end subroutine init_2d_distribution
+end subroutine init_2d_distribution
!---------------------------------------------------------------------
- elemental function heaviside(x)
+elemental function heaviside(x)
real(kind=CUSTOM_REAL), intent(in) :: x
real(kind=CUSTOM_REAL) :: heaviside
@@ -663,12 +663,12 @@
heaviside = 0e0_CUSTOM_REAL
endif
- end function heaviside
+end function heaviside
!=====================================================================
! adds boundary term Bt into Force array for each fault.
!
- subroutine bc_dynflt_set3d_all(F,Vel,Dis)
+subroutine bc_dynflt_set3d_all(F,Vel,Dis)
real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: Vel,Dis
real(kind=CUSTOM_REAL), dimension(:,:), intent(inout) :: F
@@ -684,10 +684,10 @@
endif
enddo
- end subroutine bc_dynflt_set3d_all
+end subroutine bc_dynflt_set3d_all
!---------------------------------------------------------------------
- subroutine BC_DYNFLT_set3d_TPV16(bc,MxA,V,D,iflt)
+subroutine BC_DYNFLT_set3d_TPV16(bc,MxA,V,D,iflt)
use specfem_par, only:it,NSTEP
@@ -804,13 +804,10 @@
if ( mod(it,NSNAP) == 0) call write_dataXZ(bc%dataXZ,it,iflt)
if ( it == NSTEP) call SCEC_Write_RuptureTime(bc%dataXZ,bc%dt,NSTEP,iflt)
- end subroutine BC_DYNFLT_set3d_TPV16
+end subroutine BC_DYNFLT_set3d_TPV16
-!===============================================================
-
-
!---------------------------------------------------------------------
- subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
+subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
use specfem_par, only:it,NSTEP
@@ -825,27 +822,26 @@
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
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,:))
-! get predicted values
+ ! get predicted values
dD = get_jump(bc,D) ! dD_predictor
dV = get_jump(bc,V) ! dV_predictor
dA = get_weighted_jump(bc,MxA) ! dA_free
-! rotate to fault frame (tangent,normal)
-! component 3 is normal to the fault
+ ! rotate to fault frame (tangent,normal)
+ ! component 3 is normal to the fault
dD = rotate(bc,dD,1)
dV = rotate(bc,dV,1)
dA = rotate(bc,dA,1)
-! T_stick
- T(1,:) = bc%Z * ( dV(1,:) + half_dt*dA(1,:) )
- T(2,:) = bc%Z * ( dV(2,:) + half_dt*dA(2,:) )
- T(3,:) = bc%Z * ( dV(3,:) + half_dt*dA(3,:) )
+ ! T_stick
+ T(1,:) = bc%Z * ( dV(1,:) + half_dt*dA(1,:) )
+ T(2,:) = bc%Z * ( dV(2,:) + half_dt*dA(2,:) )
+ T(3,:) = bc%Z * ( dV(3,:) + half_dt*dA(3,:) )
!Warning : dirty particular free surface condition z = 0.
! where (bc%zcoord(:) > - SMALLVAL) T(2,:) = 0
@@ -853,28 +849,28 @@
! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) < SMALLVAL) T(2,k) = 0.e0_CUSTOM_REAL
! end do
-! add initial stress
+ ! add initial stress
T = T + bc%T0
-! Solve for normal stress (negative is compressive)
+ ! Solve for normal stress (negative is compressive)
! Opening implies free stress
- if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
+ if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
-! Update slip weakening friction:
- ! Update slip state variable
- ! WARNING: during opening the friction state variable should not evolve
+ ! Update slip weakening friction:
+ ! Update slip state variable
+ ! WARNING: during opening the friction state variable should not evolve
call swf_update_state(bc%D,dD,bc%V,bc%swf)
- ! Update friction coeficient
+ ! Update friction coeficient
bc%MU = swf_mu(bc%swf)
-! combined with time-weakening for nucleation
-! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
+ ! combined with time-weakening for nucleation
+ ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
-! Update strength
+ ! Update strength
strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL)
-! Solve for shear stress
+ ! Solve for shear stress
tnorm = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
tnorm = max(tnorm,1e0_CUSTOM_REAL)
t1 = T(1,:)/tnorm
@@ -884,25 +880,25 @@
T(1,:) = tnew * t1
T(2,:) = tnew * t2
-! Save total tractions
+ ! Save total tractions
bc%T = T
-! Subtract initial stress
+ ! Subtract initial stress
T = T - bc%T0
-! Update slip acceleration da=da_free-T/(0.5*dt*Z)
+ ! Update slip acceleration da=da_free-T/(0.5*dt*Z)
dA(1,:) = dA(1,:) - T(1,:)/(bc%Z*half_dt)
dA(2,:) = dA(2,:) - T(2,:)/(bc%Z*half_dt)
dA(3,:) = dA(3,:) - T(3,:)/(bc%Z*half_dt)
-! Update slip and slip rate, in fault frame
+ ! Update slip and slip rate, in fault frame
bc%D = dD
bc%V = dV + half_dt*dA
-! Rotate tractions back to (x,y,z) frame
+ ! Rotate tractions back to (x,y,z) frame
T = rotate(bc,T,-1)
-! Add boundary term B*T to M*a
+ ! Add boundary term B*T to M*a
MxA(1,bc%ibulk1) = MxA(1,bc%ibulk1) + bc%B*T(1,:)
MxA(2,bc%ibulk1) = MxA(2,bc%ibulk1) + bc%B*T(2,:)
MxA(3,bc%ibulk1) = MxA(3,bc%ibulk1) + bc%B*T(3,:)
@@ -911,54 +907,52 @@
MxA(2,bc%ibulk2) = MxA(2,bc%ibulk2) - bc%B*T(2,:)
MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
-
-!-- intermediate storage of outputs --
+ !-- intermediate storage of outputs --
Vnorm = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
call store_dataXZ(bc%dataXZ, strength, theta_old, bc%swf%theta, bc%swf%dc, &
Vnorm_old, Vnorm, it*bc%dt,bc%dt)
call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
-
-!-- outputs --
-! write dataT every NTOUT time step or at the end of simulation
+ !-- outputs --
+ ! write dataT every NTOUT time step or at the end of simulation
if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it)
-! write dataXZ every NSNAP time step
+ ! write dataXZ every NSNAP time step
if ( mod(it,NSNAP) == 0) call write_dataXZ(bc%dataXZ,it,iflt)
if ( it == NSTEP) call SCEC_Write_RuptureTime(bc%dataXZ,bc%dt,NSTEP,iflt)
- end subroutine BC_DYNFLT_set3d
+end subroutine BC_DYNFLT_set3d
!===============================================================
- function get_jump (bc,v) result(dv)
+function get_jump (bc,v) result(dv)
type(bc_dynflt_type), intent(in) :: bc
real(kind=CUSTOM_REAL), intent(in) :: v(:,:)
real(kind=CUSTOM_REAL) :: dv(3,bc%nglob)
-! diference between side 2 and side 1 of fault nodes. dv
- dv(1,:) = v(1,bc%ibulk2)-v(1,bc%ibulk1)
- dv(2,:) = v(2,bc%ibulk2)-v(2,bc%ibulk1)
- dv(3,:) = v(3,bc%ibulk2)-v(3,bc%ibulk1)
+ ! diference between side 2 and side 1 of fault nodes. dv
+ dv(1,:) = v(1,bc%ibulk2)-v(1,bc%ibulk1)
+ dv(2,:) = v(2,bc%ibulk2)-v(2,bc%ibulk1)
+ dv(3,:) = v(3,bc%ibulk2)-v(3,bc%ibulk1)
- end function get_jump
+end function get_jump
!---------------------------------------------------------------------
- function get_weighted_jump (bc,f) result(da)
+function get_weighted_jump (bc,f) result(da)
- type(bc_dynflt_type), intent(in) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: f(:,:)
+ type(bc_dynflt_type), intent(in) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: f(:,:)
- real(kind=CUSTOM_REAL) :: da(3,bc%nglob)
+ real(kind=CUSTOM_REAL) :: da(3,bc%nglob)
-! diference between side 2 and side 1 of fault nodes. M-1 * F
- da(1,:) = bc%invM2*f(1,bc%ibulk2)-bc%invM1*f(1,bc%ibulk1)
- da(2,:) = bc%invM2*f(2,bc%ibulk2)-bc%invM1*f(2,bc%ibulk1)
- da(3,:) = bc%invM2*f(3,bc%ibulk2)-bc%invM1*f(3,bc%ibulk1)
+ ! diference between side 2 and side 1 of fault nodes. M-1 * F
+ da(1,:) = bc%invM2*f(1,bc%ibulk2)-bc%invM1*f(1,bc%ibulk1)
+ da(2,:) = bc%invM2*f(2,bc%ibulk2)-bc%invM1*f(2,bc%ibulk1)
+ da(3,:) = bc%invM2*f(3,bc%ibulk2)-bc%invM1*f(3,bc%ibulk1)
- end function get_weighted_jump
+end function get_weighted_jump
!----------------------------------------------------------------------
- function rotate(bc,v,fb) result(vr)
+function rotate(bc,v,fb) result(vr)
type(bc_dynflt_type), intent(in) :: bc
real(kind=CUSTOM_REAL), intent(in) :: v(3,bc%nglob)
@@ -981,11 +975,11 @@
endif
- end function rotate
+end function rotate
!=====================================================================
- subroutine swf_update_state(dold,dnew,vold,f)
+subroutine swf_update_state(dold,dnew,vold,f)
real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: vold,dold,dnew
type(swf_type), intent(inout) :: f
@@ -1002,11 +996,11 @@
if (vnorm<V_HEALING) f%theta(k) = 0e0_CUSTOM_REAL
enddo
endif
- end subroutine swf_update_state
+end subroutine swf_update_state
!=====================================================================
! Friction coefficient
- function swf_mu_TPV16(f,time,nglob) result(mu)
+function swf_mu_TPV16(f,time,nglob) result(mu)
type(swf_type), intent(in) :: f
real(kind=CUSTOM_REAL) :: mu(size(f%theta))
@@ -1015,36 +1009,31 @@
integer :: i,nglob
!-- linear slip weakening:
- mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
- do i=1,nglob
- if ( (f%theta(i) >= f%Dc(i)) .or. (time >= f%T(i)) ) then
- mu(i) = f%mud(i)
- endif
- enddo
+ mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
+ do i=1,nglob
+ if ( (f%theta(i) >= f%Dc(i)) .or. (time >= f%T(i)) ) then
+ mu(i) = f%mud(i)
+ endif
+ enddo
- end function swf_mu_TPV16
-
+end function swf_mu_TPV16
-
!=====================================================================
! Friction coefficient
- function swf_mu(f) result(mu)
+function swf_mu(f) result(mu)
type(swf_type), intent(in) :: f
real(kind=CUSTOM_REAL) :: mu(size(f%theta))
!-- linear slip weakening:
-
- mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
- mu = max( mu, f%mud)
+ mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
+ mu = max( mu, f%mud)
- end function swf_mu
+end function swf_mu
-
!===============================================================
! OUTPUTS
-
- subroutine init_dataT(DataT,coord,nglob,NT,iflt)
+subroutine init_dataT(DataT,coord,nglob,NT,iflt)
! NT = total number of time steps
integer, intent(in) :: nglob,NT,iflt
@@ -1115,11 +1104,10 @@
close(IIN)
- end subroutine init_dataT
+end subroutine init_dataT
-
!---------------------------------------------------------------
- subroutine store_dataT(dataT,d,v,t,itime)
+subroutine store_dataT(dataT,d,v,t,itime)
type(dataT_type), intent(inout) :: dataT
real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
@@ -1138,11 +1126,10 @@
dataT%t3(itime,i) = t(3,k)
enddo
- end subroutine store_dataT
+end subroutine store_dataT
-
!-----------------------------------------------------------------
- subroutine write_dataT_all(nt)
+subroutine write_dataT_all(nt)
integer, intent(in) :: nt
@@ -1153,10 +1140,10 @@
call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt)
enddo
- end subroutine write_dataT_all
+end subroutine write_dataT_all
!------------------------------------------------------------------------
- subroutine SCEC_write_dataT(dataT,DT,NT)
+subroutine SCEC_write_dataT(dataT,DT,NT)
type(dataT_type), intent(in) :: dataT
real(kind=CUSTOM_REAL), intent(in) :: DT
@@ -1164,14 +1151,11 @@
integer :: i,k,IOUT
character :: NTchar*5
-
-
integer :: today(3), now(3)
call idate(today) ! today(1)=day, (2)=month, (3)=year
call itime(now) ! now(1)=hour, (2)=minute, (3)=second
-
IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
write(NTchar,1) NT
@@ -1214,11 +1198,10 @@
enddo
- end subroutine SCEC_write_dataT
+end subroutine SCEC_write_dataT
!-------------------------------------------------------------------------------------------------
-
- subroutine SCEC_Write_RuptureTime(dataXZ,DT,NT,iflt)
+subroutine SCEC_Write_RuptureTime(dataXZ,DT,NT,iflt)
type(dataXZ_type), intent(in) :: dataXZ
real(kind=CUSTOM_REAL), intent(in) :: DT
@@ -1226,7 +1209,6 @@
integer :: i,IOUT
character(len=70) :: filename
-
integer*4 today(3), now(3)
call idate(today) ! today(1)=day, (2)=month, (3)=year
@@ -1257,14 +1239,11 @@
1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
+end subroutine SCEC_Write_RuptureTime
-
- end subroutine SCEC_Write_RuptureTime
-
!-------------------------------------------------------------------------------------------------
+subroutine init_dataXZ(DataXZ,bc,nglob)
- subroutine init_dataXZ(DataXZ,bc,nglob)
-
type(dataXZ_type), intent(inout) :: DataXZ
type(bc_dynflt_type) :: bc
integer, intent(in) :: nglob
@@ -1288,9 +1267,8 @@
DataXZ%tRUP = 0e0_CUSTOM_REAL
DataXZ%tPZ = 0e0_CUSTOM_REAL
+end subroutine init_dataXZ
- end subroutine init_dataXZ
-
!---------------------------------------------------------------
subroutine store_dataXZ(dataXZ,stg,dold,dnew,dc,vold,vnew,time,dt)
@@ -1326,18 +1304,16 @@
! note: the other arrays in dataXZ are pointers to arrays in bc
! they do not need to be updated here
- end subroutine store_dataXZ
+end subroutine store_dataXZ
!---------------------------------------------------------------
- subroutine write_dataXZ(dataXZ,itime,iflt)
+subroutine write_dataXZ(dataXZ,itime,iflt)
-
type(dataXZ_type), intent(in) :: dataXZ
integer, intent(in) :: itime,iflt
character(len=70) :: filename
-
write(filename,"('OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
! 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
@@ -1363,7 +1339,7 @@
write(IOUT) dataXZ%tPZ
close(IOUT)
- end subroutine write_dataXZ
+end subroutine write_dataXZ
end module fault_solver
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-02-25 01:39:11 UTC (rev 19681)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2012-02-25 16:39:12 UTC (rev 19682)
@@ -337,7 +337,7 @@
! (nint : fortran round (nearest whole number) ,
! if nint(a)=0.5 then "a" get upper bound )
-! Loading the next slipt_rate one ahead it.
+! Loading the next slip_rate one ahead it.
! This is done in case bc%kin_dt
! if (it_kin == it) it_kin=it_kin+1 !
More information about the CIG-COMMITS
mailing list