[cig-commits] r20931 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
ampuero at geodynamics.org
ampuero at geodynamics.org
Thu Oct 25 17:16:57 PDT 2012
Author: ampuero
Date: 2012-10-25 17:16:57 -0700 (Thu, 25 Oct 2012)
New Revision: 20931
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90
Log:
some clean up in solver_common
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-10-25 19:35:03 UTC (rev 20930)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-10-26 00:16:57 UTC (rev 20931)
@@ -319,7 +319,7 @@
enddo
-subroutine TPV16_init
+end subroutine TPV16_init
end subroutine init_one_fault
@@ -666,7 +666,7 @@
mu = swf_mu(f)
-subroutine swf_init
+end subroutine swf_init
!---------------------------------------------------------------------
subroutine swf_update_state(dold,dnew,vold,f)
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90 2012-10-25 19:35:03 UTC (rev 20930)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90 2012-10-26 00:16:57 UTC (rev 20931)
@@ -42,11 +42,9 @@
real(kind=CUSTOM_REAL) :: tmp_vec(3,NGLOB_AB)
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: jacobian2Dw
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: normal
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: nxyz
integer, dimension(:,:), allocatable :: ibool1
- real(kind=CUSTOM_REAL) :: norm
- integer :: n1,n2,n3
integer :: ij,k,e
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: nx,ny,nz
read(IIN_BIN) bc%nspec,bc%nglob
if (.NOT.PARALLEL_FAULT .and. bc%nspec==0) return
@@ -77,16 +75,12 @@
bc%dt = dt
bc%B = 0e0_CUSTOM_REAL
- allocate(nx(bc%nglob),ny(bc%nglob),nz(bc%nglob))
- nx = 0e0_CUSTOM_REAL
- ny = 0e0_CUSTOM_REAL
- nz = 0e0_CUSTOM_REAL
+ allocate(nxyz(3,bc%nglob))
+ nxyz = 0e0_CUSTOM_REAL
do e=1,bc%nspec
do ij = 1,NGLLSQUARE
k = ibool1(ij,e)
- nx(k) = nx(k) + normal(1,ij,e)
- ny(k) = ny(k) + normal(2,ij,e)
- nz(k) = nz(k) + normal(3,ij,e)
+ nxyz(:,k) = nxyz(:,k) + normal(:,ij,e)
bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
enddo
enddo
@@ -104,34 +98,20 @@
if (bc%nspec>0) bc%B = tmp_vec(1,bc%ibulk1)
tmp_vec = 0._CUSTOM_REAL
- if (bc%nspec>0) then
- tmp_vec(1,bc%ibulk1) = nx
- tmp_vec(2,bc%ibulk1) = ny
- tmp_vec(3,bc%ibulk1) = nz
- endif
+ if (bc%nspec>0) tmp_vec(:,bc%ibulk1) = nxyz
! assembles with other MPI processes
call assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,tmp_vec, &
num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
my_neighbours_ext_mesh)
- if (bc%nspec>0) then
- nx = tmp_vec(1,bc%ibulk1)
- ny = tmp_vec(2,bc%ibulk1)
- nz = tmp_vec(3,bc%ibulk1)
- endif
+ if (bc%nspec>0) nxyz = tmp_vec(:,bc%ibulk1)
endif
if (bc%nspec>0) then
- do k=1,bc%nglob
- norm = sqrt( nx(k)*nx(k) + ny(k)*ny(k) + nz(k)*nz(k) )
- nx(k) = nx(k) / norm
- ny(k) = ny(k) / norm
- nz(k) = nz(k) / norm
- enddo
+ call normalize_3d_vector(nxyz)
+ call compute_R(bc%R,bc%nglob,nxyz)
- call compute_R(bc%R,bc%nglob,nx,ny,nz)
-
! Needed in dA_Free = -K2*d2/M2 + K1*d1/M1
bc%invM1 = Minv(bc%ibulk1)
bc%invM2 = Minv(bc%ibulk2)
@@ -150,38 +130,52 @@
end subroutine initialize_fault
!---------------------------------------------------------------------
-subroutine compute_R(R,nglob,nx,ny,nz)
+subroutine normalize_3d_vector(v)
+ real(kind=CUSTOM_REAL), intent(inout) :: v(:,:)
+
+ real(kind=CUSTOM_REAL) :: norm
+ integer :: k
+
+ ! assume size(v) = [3,N]
+ do k=1,size(v,2)
+ norm = sqrt( v(1,k)*v(1,k) + v(2,k)*v(2,k) + v(3,k)*v(3,k) )
+ v(:,k) = v(:,k) / norm
+ enddo
+
+end subroutine normalize_3d_vector
+
+!---------------------------------------------------------------------
+! Percy: define fault directions according to SCEC conventions
+! Fault coordinates (s,d,n) = (1,2,3)
+! s = strike , d = dip , n = normal
+! 1 = strike , 2 = dip , 3 = normal
+! with dip pointing downwards
+!
+subroutine compute_R(R,nglob,n)
+
integer :: nglob
real(kind=CUSTOM_REAL), intent(out) :: R(3,3,nglob)
- real(kind=CUSTOM_REAL), dimension(nglob), intent(in) :: nx,ny,nz
+ real(kind=CUSTOM_REAL), intent(in) :: n(3,nglob)
- real(kind=CUSTOM_REAL), dimension(nglob) :: sx,sy,sz,dx,dy,dz,norm
+ real(kind=CUSTOM_REAL), dimension(3,nglob) :: s,d
- ! Percy , defining fault directions (in concordance with SCEC conventions) .
- ! fault coordinates (s,d,n) = (1,2,3)
- ! s = strike , d = dip , n = n.
- ! 1 = strike , 2 = dip , 3 = n.
- norm = sqrt(nx*nx+ny*ny)
- sx = ny/norm
- sy = -nx/norm
- sz = 0.e0_CUSTOM_REAL
+ s(1,:) = n(2,:) ! sx = ny
+ s(2,:) = -n(1,:) ! sy =-nx
+ s(3,:) = 0.e0_CUSTOM_REAL
+ call normalize_3d_vector(s)
- norm = sqrt(sy*sy*nz*nz+sx*sx*nz*nz+(sy*nx-ny*sx)*(nx*sy-ny*sx))
- dx = -sy*nz/norm
- dy = sx*nz/norm
- dz = (sy*nx-ny*sx)/norm
- !Percy, dz is always dipwards = -1/norm , because (nx*sy-ny*sx)= - 1
+ d(1,:) = -s(2,:)*n(3,:) ! dx = -sy*nz
+ d(2,:) = s(1,:)*n(3,:) ! dy = sx*nz
+ d(3,:) = s(2,:)*n(1,:) - s(1,:)*n(2,:) ! dz = sy*nx-ny*sx
+ call normalize_3d_vector(d)
+ ! dz is always dipwards (negative), because
+ ! (nx*sy-ny*sx) = -(nx^2+ny^2)/sqrt(nx^2+ny^2)
+ ! = -sqrt(nx^2+ny^2) < 0
- R(1,1,:)=sx
- R(1,2,:)=sy
- R(1,3,:)=sz
- R(2,1,:)=dx
- R(2,2,:)=dy
- R(2,3,:)=dz
- R(3,1,:)=nx
- R(3,2,:)=ny
- R(3,3,:)=nz
+ R(1,:,:) = s
+ R(2,:,:) = d
+ R(3,:,:) = n
end subroutine compute_R
More information about the CIG-COMMITS
mailing list