[cig-commits] r20902 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Wed Oct 24 07:09:25 PDT 2012
Author: surendra
Date: 2012-10-24 07:09:25 -0700 (Wed, 24 Oct 2012)
New Revision: 20902
Modified:
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:
Parallelized kinematic fault solver
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-24 13:48:12 UTC (rev 20901)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-10-24 14:09:25 UTC (rev 20902)
@@ -617,15 +617,6 @@
! call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
call init_dataXZ(bc%dataXZ,bc)
-!JPA for debugging: test write_dataXZ by printing a snapshot at it=0 with processor id in place of d1
- if (PARALLEL_FAULT) then
- if(bc%nglob > 0) bc%dataXZ%d1 = myrank
- call gather_dataXZ(bc)
- if (myrank==0) call write_dataXZ(bc%dataXZ_all,0,iflt)
- if(bc%nglob > 0) bc%dataXZ%d1 = 0.e0_CUSTOM_REAL
- endif
-
-
end subroutine init_one_fault
!=====================================================================
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-10-24 13:48:12 UTC (rev 20901)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2012-10-24 14:09:25 UTC (rev 20902)
@@ -139,128 +139,128 @@
NAMELIST / KINPAR / kindt
read(IIN_BIN) bc%nspec,bc%nglob
- if (bc%nspec==0) return
+ if (.NOT.PARALLEL_FAULT .and. bc%nspec==0) return
+ if (bc%nspec>0) then
- allocate( bc%ibulk1(bc%nglob) )
- allocate( bc%ibulk2(bc%nglob) )
- allocate( ibool1(NGLLSQUARE,bc%nspec) )
- allocate(normal(NDIM,NGLLSQUARE,bc%nspec))
- allocate(jacobian2Dw(NGLLSQUARE,bc%nspec))
- allocate(bc%coord(3,bc%nglob))
+ allocate( bc%ibulk1(bc%nglob) )
+ allocate( bc%ibulk2(bc%nglob) )
+ allocate( ibool1(NGLLSQUARE,bc%nspec) )
+ allocate(normal(NDIM,NGLLSQUARE,bc%nspec))
+ allocate(jacobian2Dw(NGLLSQUARE,bc%nspec))
+ allocate(bc%coord(3,bc%nglob))
- read(IIN_BIN) ibool1
- read(IIN_BIN) jacobian2Dw
- read(IIN_BIN) normal
- read(IIN_BIN) bc%ibulk1
- read(IIN_BIN) bc%ibulk2
- read(IIN_BIN) bc%coord(1,:)
- read(IIN_BIN) bc%coord(2,:)
- read(IIN_BIN) bc%coord(3,:)
- bc%dt = dt
+ read(IIN_BIN) ibool1
+ read(IIN_BIN) jacobian2Dw
+ read(IIN_BIN) normal
+ read(IIN_BIN) bc%ibulk1
+ read(IIN_BIN) bc%ibulk2
+ read(IIN_BIN) bc%coord(1,:)
+ read(IIN_BIN) bc%coord(2,:)
+ read(IIN_BIN) bc%coord(3,:)
+ bc%dt = dt
- allocate( bc%B(bc%nglob) )
- 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
- 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)
- bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
+ allocate( bc%B(bc%nglob) )
+ 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
+ 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)
+ bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
+ enddo
enddo
- enddo
- ! TO DO: assemble B and n across processors
- 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
- allocate( bc%R(3,3,bc%nglob) )
- call compute_R(bc%R,bc%nglob,nx,ny,nz)
- deallocate(nx,ny,nz)
- ! inverse M needed in dA_Free = -K2*d2/M2 + K1*d1/M1
- allocate(bc%invM1(bc%nglob))
- allocate(bc%invM2(bc%nglob))
- bc%invM1 = Minv(bc%ibulk1)
- bc%invM2 = Minv(bc%ibulk2)
+ endif
- ! Fault impedance, Z in : Trac=T_Stick-Z*dV
- ! Z = 1/( B1/M1 + B2/M2 ) / (0.5*dt)
- ! T_Stick = Z*Vfree traction as if the fault was stuck (no displ discontinuity)
- ! NOTE: same Bi on both sides, see note above
- allocate(bc%Z(bc%nglob))
- bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
- ! WARNING: In non-split nodes at fault edges Minv is assembled across the fault
- ! hence invM1+invM2=2/(M1+M2) instead of 1/M1+1/M2
- ! For instance, in a symmetric mesh (M1=M2) Z is twice its intended value
+ if (PARALLEL_FAULT) then
- allocate(bc%T(3,bc%nglob))
- allocate(bc%D(3,bc%nglob))
- allocate(bc%V(3,bc%nglob))
- bc%T = 0e0_CUSTOM_REAL
- bc%D = 0e0_CUSTOM_REAL
- bc%V = 0e0_CUSTOM_REAL
-
- ! time interval between two loaded slip rates
- read(IIN_PAR,nml=KINPAR)
- bc%kin_dt = kindt
-
- bc%kin_it=0
- ! Always have in memory the slip-rate model at two times, t1 and t2,
- ! spatially interpolated in the spectral element grid
- allocate(bc%v_kin_t1(2,bc%nglob))
- allocate(bc%v_kin_t2(2,bc%nglob))
- bc%v_kin_t1 = 0e0_CUSTOM_REAL
- bc%v_kin_t2 = 0e0_CUSTOM_REAL
+ accel=0._CUSTOM_REAL
+ if (bc%nspec>0) accel(1,bc%ibulk1) = bc%B
+ ! assembles with other MPI processes
+ call assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,accel, &
+ 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) bc%B = accel(1,bc%ibulk1)
+
+ accel=0._CUSTOM_REAL
+ if (bc%nspec>0) then
+ accel(1,bc%ibulk1) = nx
+ accel(2,bc%ibulk1) = ny
+ accel(3,bc%ibulk1) = nz
+ endif
+ ! assembles with other MPI processes
+ call assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,accel, &
+ 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 = accel(1,bc%ibulk1)
+ ny = accel(2,bc%ibulk1)
+ nz = accel(3,bc%ibulk1)
+ endif
- call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
- call init_dataXZ(bc%dataXZ,bc%nglob)
+ endif
-end subroutine init_one_fault
+ 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
+ allocate( bc%R(3,3,bc%nglob) )
+ call compute_R(bc%R,bc%nglob,nx,ny,nz)
+ deallocate(nx,ny,nz)
-!---------------------------------------------------------------------
-subroutine compute_R(R,nglob,nx,ny,nz)
+ ! inverse M needed in dA_Free = -K2*d2/M2 + K1*d1/M1
+ allocate(bc%invM1(bc%nglob))
+ allocate(bc%invM2(bc%nglob))
+ bc%invM1 = Minv(bc%ibulk1)
+ bc%invM2 = Minv(bc%ibulk2)
- integer :: nglob
- real(kind=CUSTOM_REAL), intent(out) :: R(3,3,nglob)
- real(kind=CUSTOM_REAL), dimension(nglob), intent(in) :: nx,ny,nz
+ ! Fault impedance, Z in : Trac=T_Stick-Z*dV
+ ! Z = 1/( B1/M1 + B2/M2 ) / (0.5*dt)
+ ! T_Stick = Z*Vfree traction as if the fault was stuck (no displ discontinuity)
+ ! NOTE: same Bi on both sides, see note above
+ allocate(bc%Z(bc%nglob))
+ bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
+ ! WARNING: In non-split nodes at fault edges Minv is assembled across the fault
+ ! hence invM1+invM2=2/(M1+M2) instead of 1/M1+1/M2
+ ! For instance, in a symmetric mesh (M1=M2) Z is twice its intended value
- real(kind=CUSTOM_REAL), dimension(nglob) :: sx,sy,sz,dx,dy,dz,norm
+ allocate(bc%T(3,bc%nglob))
+ allocate(bc%D(3,bc%nglob))
+ allocate(bc%V(3,bc%nglob))
+ bc%T = 0e0_CUSTOM_REAL
+ bc%D = 0e0_CUSTOM_REAL
+ bc%V = 0e0_CUSTOM_REAL
+
+ ! time interval between two loaded slip rates
+ read(IIN_PAR,nml=KINPAR)
+ bc%kin_dt = kindt
+
+ bc%kin_it=0
+ ! Always have in memory the slip-rate model at two times, t1 and t2,
+ ! spatially interpolated in the spectral element grid
+ allocate(bc%v_kin_t1(2,bc%nglob))
+ allocate(bc%v_kin_t2(2,bc%nglob))
+ bc%v_kin_t1 = 0e0_CUSTOM_REAL
+ bc%v_kin_t2 = 0e0_CUSTOM_REAL
- ! Percy: defining fault directions following 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
+ call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
+ call init_dataXZ(bc%dataXZ,bc%nglob)
- 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
+ endif
- 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
+end subroutine init_one_fault
-end subroutine compute_R
-
!=====================================================================
! adds boundary term Bt to Force array for each fault.
!
@@ -299,103 +299,105 @@
real(kind=CUSTOM_REAL) :: t1,t2
real(kind=CUSTOM_REAL) :: half_dt,time
- half_dt = 0.5e0_CUSTOM_REAL*bc%dt
+ if (bc%nspec > 0) then !Surendra : for parallel faults
- ! 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
+ half_dt = 0.5e0_CUSTOM_REAL*bc%dt
- ! 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)
+ ! 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
- ! Time marching
- time = it*bc%dt
- ! Slip_rate step "it_kin"
- it_kin = bc%kin_it*nint(bc%kin_dt/bc%dt)
- ! (nint : fortran round (nearest whole number) ,
- ! if nint(a)=0.5 then "a" get upper bound )
+ ! 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)
- ! 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 !
+ ! Time marching
+ time = it*bc%dt
+ ! Slip_rate step "it_kin"
+ it_kin = bc%kin_it*nint(bc%kin_dt/bc%dt)
+ ! (nint : fortran round (nearest whole number) ,
+ ! if nint(a)=0.5 then "a" get upper bound )
+ ! 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 !
- !NOTE : it and it_kin is being used due to integers are exact numbers.
- if (it > it_kin) then
- print*, 'it :'
- print*, it
- print*, 'it_kin'
- print*, it_kin
+ !NOTE : it and it_kin is being used due to integers are exact numbers.
+ if (it > it_kin) then
- bc%kin_it = bc%kin_it +1
- bc%v_kin_t1 = bc%v_kin_t2
- print*, 'loading v_kin_t2'
- !Temporal : just for snapshots file names kin_dt=0.1 , dt=0.0001
- !snapshot(100=itime).. : itime=kin_it*(kin_dt/dt)
- itime = bc%kin_it*nint(bc%kin_dt/bc%dt)
- call load_vslip_snapshots(bc%dataXZ,itime,bc%nglob,iflt)
-! loading slip rates
- bc%v_kin_t2(1,:)=bc%dataXZ%v1
- bc%v_kin_t2(2,:)=bc%dataXZ%v2
+ print*, 'it :', it
+ print*, 'it_kin :', it_kin
+
+ bc%kin_it = bc%kin_it +1
+ bc%v_kin_t1 = bc%v_kin_t2
+ print*, 'loading v_kin_t2'
+ !Temporal : just for snapshots file names kin_dt=0.1 , dt=0.0001
+ !snapshot(100=itime).. : itime=kin_it*(kin_dt/dt)
+ itime = bc%kin_it*nint(bc%kin_dt/bc%dt)
+ call load_vslip_snapshots(bc%dataXZ,itime,bc%nglob,iflt)
+! loading slip rates
+ bc%v_kin_t2(1,:)=bc%dataXZ%v1
+ bc%v_kin_t2(2,:)=bc%dataXZ%v2
- !linear interpolation in time between t1 and t2
- !REMARK , bc%kin_dt is the delta "t" between two snapshots.
+ !linear interpolation in time between t1 and t2
+ !REMARK , bc%kin_dt is the delta "t" between two snapshots.
- endif
+ endif
- t1 = (bc%kin_it-1) * bc%kin_dt
- t2 = bc%kin_it * bc%kin_dt
+ t1 = (bc%kin_it-1) * bc%kin_dt
+ t2 = bc%kin_it * bc%kin_dt
- ! Kinematic velocity_rate
- ! bc%V : Imposed a priori and read from slip rate snapshots (from time reversal)
- ! Linear interpolation between consecutive kinematic time steps.
- ! V will be given at 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
+ ! Kinematic velocity_rate
+ ! bc%V : Imposed a priori and read from slip rate snapshots (from time reversal)
+ ! Linear interpolation between consecutive kinematic time steps.
+ ! V will be given at 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,:)
+ !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) , V known apriori as input.
- ! CONVENTION : T(ibulk1)=T=-T(ibulk2)
- 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,:) )
+ ! T = Z*( dV_free - V) , V known apriori as input.
+ ! CONVENTION : T(ibulk1)=T=-T(ibulk2)
+ 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
+ ! Save tractions
+ bc%T = T
- ! Update slip in fault frame
- bc%D = dD
+ ! Update slip in fault frame
+ bc%D = dD
- ! Rotate tractions back to (x,y,z) frame
- T = rotate(bc,T,-1)
+ ! Rotate tractions back to (x,y,z) frame
+ T = rotate(bc,T,-1)
- ! 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,:)
+ ! 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,:)
- MxA(1,bc%ibulk2) = MxA(1,bc%ibulk2) - bc%B*T(1,:)
- MxA(2,bc%ibulk2) = MxA(2,bc%ibulk2) - bc%B*T(2,:)
- MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
+ MxA(1,bc%ibulk2) = MxA(1,bc%ibulk2) - bc%B*T(1,:)
+ 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 --
- call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
+ !-- intermediate storage of outputs --
+ 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
- if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it)
- ! write dataXZ every NSNAP time steps
- ! if ( mod(it,NSNAP) == 0) call write_dataXZ(bc,it,iflt)
+ !-- OUTPUTS --
+ ! write dataT every NTOUT time steps 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 steps
+ ! if ( mod(it,NSNAP) == 0) call write_dataXZ(bc,it,iflt)
+ endif
+
end subroutine BC_KINFLT_set_single
!===============================================================
More information about the CIG-COMMITS
mailing list