Skip to content

Commit

Permalink
Merge pull request #1693 from Huihuiweng/devel
Browse files Browse the repository at this point in the history 
Refine TWF by adding independent mus and mud and a special twf for TPV22; Update the corresponding user manual
  • Loading branch information
@danielpeter
danielpeter authored Apr 25, 2024
2 parents 7a1c764 + 8cc2623 commit ce914f7
Show file tree
Hide file tree
Showing 4 changed files with 74 additions and 29 deletions.
13 changes: 10 additions & 3 deletions doc/USER_MANUAL/06_fault_sources.tex
View file
Original file line number Diff line number Diff line change
Expand Up @@ -246,8 +246,8 @@ \section{Input Files}
followed by (n1+n2+n3) \&\textbf{DIST2D} blocks\newline


\&\textbf{SWF} mus, mud, dc {[}, nmus, nmud, ndc{]} (weakening\_kind=1 for linear (default); weakening\_kind=2 for exponential) /\newline
\&\textbf{TWF} nuc\_x, nuc\_y, nuc\_z, nuc\_r, nuc\_t0, nuc\_v /\newline
\&\textbf{SWF} mus, mud, dc {[}, nmus, nmud, ndc{]} (weakening\_kind=1 for linear (default); weakening\_kind=2 for exponential; weakening\_kind=3 for power-law with p) /\newline
\&\textbf{TWF} nuc\_x, nuc\_y, nuc\_z, nuc\_r, nuc\_t0, nuc\_v, mus, mud, kind /\newline
\&\textbf{RSF} V0,f0,a,b,L,V\_init,theta\_init,C,StateLaw {[} nV0,nf0,na,nb,nL,nV\_init,ntheta\_init,nC {]} /\newline


Expand Down Expand Up @@ -315,8 +315,12 @@ \section{Input Files}
\item [{ndc}] = number of heterogeneous items for critical slip-weakening
distance {[}default is 0{]}
\item [{nC}] = number of heterogeneous items for cohesion
\item [{p}] = power-law coefficient (Chambon et al., 2006)
{[}default is 1{]}
\item [{weakening\_kind}] = 1 for linear slip-weakening law (default); 2 for exponential slip-weakening law:
$$\mu = \mu_d + (\mu_s-\mu_d) exp(-\frac{u}{d_c}) $$
$$\mu = \mu_d + (\mu_s-\mu_d) exp(-\frac{u}{d_c}) $$;
3 for power-law slip-weakening law:
$$\mu = \mu_d + (\mu_s-\mu_d) (1+\frac{u}{p d_c})^{-p} $$
\end{description}

\item [{\&\textbf{TWF}}] input block sets the time-weakening friction parameters
Expand All @@ -328,6 +332,9 @@ \section{Input Files}
\item [{nuc\_r}] = the radius of time-weakening nucleation (in m)
\item [{nuc\_t0}] = the cohesive time. The friction linearly decreases from static friction to dynamic friction during this time (in s)
\item [{nuc\_v}] = the time-weakening nucleation speed (in m/s)
\item [{mus}] = constant static friction coefficient for TWF
\item [{mud}] = constant dynamic friction coefficient for TWF
\item [{kind}] = 1 for regular TWF (default); 2 for TWF in TPV22* benchmarks
\end{description}

\item [{\&\textbf{RSF}}] input block sets the rate and state friction parameters
Expand Down
Binary file added my_mpi.mod
View file
Binary file not shown.
2 changes: 1 addition & 1 deletion src/specfem3D/fault_solver_common.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ module fault_solver_common
end type swf_type

type twf_type
real(kind=CUSTOM_REAL) :: nuc_x, nuc_y, nuc_z, nuc_r, nuc_t0, nuc_v
real(kind=CUSTOM_REAL) :: nuc_x,nuc_y,nuc_z,nuc_r,nuc_t0,nuc_v,mus,mud,kind
end type twf_type


Expand Down
88 changes: 63 additions & 25 deletions src/specfem3D/fault_solver_dynamic.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -996,6 +996,7 @@ subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
real(kind=CUSTOM_REAL) :: deltat,half_dt,timeval
integer :: i,ipoin,iglob,ipass
real(kind=CUSTOM_REAL) :: nuc_x, nuc_y, nuc_z, nuc_r, nuc_t0, nuc_v, dist, tw_r, coh_size
real(kind=CUSTOM_REAL) :: temp_T

! note: this implementation follows the description in:
! - rate and state friction:
Expand Down Expand Up @@ -1108,18 +1109,32 @@ subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
nuc_r = bc%twf%nuc_r
nuc_t0 = bc%twf%nuc_t0
nuc_v = bc%twf%nuc_v
do i = 1,bc%nglob
dist = ((bc%coord(1,i)-nuc_x)**2 + (bc%coord(2,i)-nuc_y)**2 + (bc%coord(3,i)-nuc_z)**2)**0.5
if (dist <= nuc_r) then
tw_r = timeval * nuc_v
coh_size = nuc_t0 * nuc_v
if (dist <= tw_r - coh_size) then
bc%mu(i) = min(bc%mu(i), bc%swf%mud(i))
else if (dist > tw_r - coh_size .and. dist <= tw_r ) then
bc%mu(i) = min(bc%mu(i), bc%swf%mud(i) + (dist-(tw_r-coh_size))/coh_size*(bc%swf%mus(i)-bc%swf%mud(i)))
if(bc%twf%kind==1) then
do i = 1,bc%nglob
dist = ((bc%coord(1,i)-nuc_x)**2 + (bc%coord(2,i)-nuc_y)**2 + (bc%coord(3,i)-nuc_z)**2)**0.5
if (dist < nuc_r) then
tw_r = timeval * nuc_v
coh_size = nuc_t0 * nuc_v
if (dist <= tw_r - coh_size) then
bc%mu(i) = min(bc%mu(i), bc%twf%mud)
else if (dist > tw_r - coh_size .and. dist <= tw_r ) then
bc%mu(i) = min(bc%mu(i), bc%twf%mud + (dist-(tw_r-coh_size))/coh_size*(bc%twf%mus-bc%twf%mud))
endif
endif
endif
enddo
enddo
elseif(bc%twf%kind==2) then
do i = 1,bc%nglob
dist = ((bc%coord(1,i)-nuc_x)**2 + (bc%coord(2,i)-nuc_y)**2 + (bc%coord(3,i)-nuc_z)**2)**0.5
if (dist < nuc_r) then
temp_T = dist/0.7/3464 + 0.081*nuc_r/0.7/3464*(1/(1-(dist/nuc_r)**2)-1)
if (timeval >= temp_T+nuc_t0) then
bc%mu(i) = min(bc%mu(i), bc%twf%mud)
else if (timeval>=temp_T .and. timeval<temp_T+nuc_t0) then
bc%mu(i) = min(bc%mu(i), bc%twf%mud + (temp_T+nuc_t0-timeval)/nuc_t0*(bc%twf%mus-bc%twf%mud))
endif
endif
enddo
endif
endif

! TPV16 benchmark
Expand Down Expand Up @@ -1299,6 +1314,7 @@ subroutine BC_DYNFLT_set3d_lts(bc,MxA,V,D,iflt,p_value,ilevel)
real(kind=CUSTOM_REAL) :: deltat,half_dt,timeval
integer :: i,ipoin,iglob,ipass
real(kind=CUSTOM_REAL) :: nuc_x, nuc_y, nuc_z, nuc_r, nuc_t0, nuc_v, dist, tw_r, coh_size
real(kind=CUSTOM_REAL) :: temp_T
! LTS
integer :: inum,num_p_local,iilevel
! only update p-level nodes (when called for each ilevel separately) set this to .true., otherwise
Expand Down Expand Up @@ -1510,18 +1526,32 @@ subroutine BC_DYNFLT_set3d_lts(bc,MxA,V,D,iflt,p_value,ilevel)
nuc_r = bc%twf%nuc_r
nuc_t0 = bc%twf%nuc_t0
nuc_v = bc%twf%nuc_v
do i = 1,bc%nglob
dist = ((bc%coord(1,i)-nuc_x)**2 + (bc%coord(2,i)-nuc_y)**2 + (bc%coord(3,i)-nuc_z)**2)**0.5
if (dist <= nuc_r) then
tw_r = timeval * nuc_v
coh_size = nuc_t0 * nuc_v
if (dist <= tw_r - coh_size) then
bc%mu(i) = min(bc%mu(i), bc%swf%mud(i))
else if (dist > tw_r - coh_size .and. dist <= tw_r ) then
bc%mu(i) = min(bc%mu(i), bc%swf%mud(i) + (dist-(tw_r-coh_size))/coh_size*(bc%swf%mus(i)-bc%swf%mud(i)))
if(bc%twf%kind==1) then
do i = 1,bc%nglob
dist = ((bc%coord(1,i)-nuc_x)**2 + (bc%coord(2,i)-nuc_y)**2 + (bc%coord(3,i)-nuc_z)**2)**0.5
if (dist < nuc_r) then
tw_r = timeval * nuc_v
coh_size = nuc_t0 * nuc_v
if (dist <= tw_r - coh_size) then
bc%mu(i) = min(bc%mu(i), bc%twf%mud)
else if (dist > tw_r - coh_size .and. dist <= tw_r ) then
bc%mu(i) = min(bc%mu(i), bc%twf%mud + (dist-(tw_r-coh_size))/coh_size*(bc%twf%mus-bc%twf%mud))
endif
endif
endif
enddo
enddo
elseif(bc%twf%kind==2) then
do i = 1,bc%nglob
dist = ((bc%coord(1,i)-nuc_x)**2 + (bc%coord(2,i)-nuc_y)**2 + (bc%coord(3,i)-nuc_z)**2)**0.5
if (dist < nuc_r) then
temp_T = dist/0.7/3464 + 0.081*nuc_r/0.7/3464*(1/(1-(dist/nuc_r)**2)-1)
if (timeval >= temp_T+nuc_t0) then
bc%mu(i) = min(bc%mu(i), bc%twf%mud)
else if (timeval>=temp_T .and. timeval<temp_T+nuc_t0) then
bc%mu(i) = min(bc%mu(i), bc%twf%mud + (temp_T+nuc_t0-timeval)/nuc_t0*(bc%twf%mus-bc%twf%mud))
endif
endif
enddo
endif
endif

! TPV16 benchmark
Expand Down Expand Up @@ -1860,7 +1890,7 @@ subroutine swf_init(f,mu,coord,IIN_PAR)
mus = 0.6_CUSTOM_REAL
mud = 0.1_CUSTOM_REAL
dc = 1.0_CUSTOM_REAL
p = 0.0_CUSTOM_REAL
p = 1.0_CUSTOM_REAL
C = 0.0_CUSTOM_REAL
T = HUGEVAL

Expand Down Expand Up @@ -1906,8 +1936,8 @@ subroutine twf_init(f,IIN_PAR)

integer :: ier

real(kind=CUSTOM_REAL) :: nuc_x, nuc_y, nuc_z, nuc_r, nuc_t0, nuc_v
NAMELIST / TWF / nuc_x, nuc_y, nuc_z, nuc_r,nuc_t0,nuc_v
real(kind=CUSTOM_REAL) :: nuc_x, nuc_y, nuc_z, nuc_r, nuc_t0, nuc_v, mus, mud, kind
NAMELIST / TWF / nuc_x, nuc_y, nuc_z, nuc_r,nuc_t0,nuc_v,mus,mud,kind

nuc_x = 0.0_CUSTOM_REAL
nuc_y = 0.0_CUSTOM_REAL
Expand All @@ -1917,6 +1947,10 @@ subroutine twf_init(f,IIN_PAR)
nuc_t0 = 0.0_CUSTOM_REAL
nuc_v = 0.0_CUSTOM_REAL

mus = 0.0_CUSTOM_REAL
mud = 0.0_CUSTOM_REAL
kind = 1

read(IIN_PAR, nml=TWF,iostat=ier)
if (ier /= 0) write(*,*) 'TWF not found in Par_file_faults.'

Expand All @@ -1928,6 +1962,10 @@ subroutine twf_init(f,IIN_PAR)
f%nuc_t0 = nuc_t0
f%nuc_v = nuc_v

f%mus = mus
f%mud = mud
f%kind = kind

end subroutine twf_init


Expand Down

0 comments on commit ce914f7

Please sign in to comment.