[cig-commits] r21135 - in seismo/3D/SPECFEM3D/trunk/src: decompose_mesh specfem3D
surendra at geodynamics.org
surendra at geodynamics.org
Tue Dec 11 18:03:01 PST 2012
Author: surendra
Date: 2012-12-11 18:03:01 -0800 (Tue, 11 Dec 2012)
New Revision: 21135
Modified:
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/fault_scotch.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_common.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_dynamic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_kinematic.f90
Log:
Got rid of warning messages.
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/fault_scotch.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/fault_scotch.f90 2012-12-12 01:57:36 UTC (rev 21134)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/fault_scotch.f90 2012-12-12 02:03:01 UTC (rev 21135)
@@ -209,24 +209,22 @@
! Lexicographic reordering of fault elements (based on their centroid)
! to make sure both sides are ordered in the same way
! and hence elements facing each other have the same index
- subroutine reorder_fault_elements(nodes_coords,nnodes)
+ subroutine reorder_fault_elements(nodes_coords)
- integer, intent(in) :: nnodes
double precision,dimension(3,nnodes), intent(in) :: nodes_coords
integer :: i
do i=1,size(faults)
- call reorder_fault_elements_single(faults(i),nodes_coords,nnodes)
+ call reorder_fault_elements_single(faults(i),nodes_coords)
enddo
end subroutine reorder_fault_elements
! ---------------------------------------------------------------------------------------------------
- subroutine reorder_fault_elements_single(f,nodes_coords,nnodes)
+ subroutine reorder_fault_elements_single(f,nodes_coords)
type(fault_type), intent(inout) :: f
- integer, intent(in) :: nnodes
double precision, dimension(3,nnodes), intent(in) :: nodes_coords
double precision, dimension(3,f%nspec) :: xyz_c
@@ -243,7 +241,7 @@
enddo
xyz_c = xyz_c / 4d0
! reorder
- call lex_order(xyz_c,new_index_list,nnodes,f%nspec)
+ call lex_order(xyz_c,new_index_list,f%nspec)
f%ispec1 = f%ispec1(new_index_list)
f%inodes1 = f%inodes1(:,new_index_list)
@@ -256,16 +254,16 @@
enddo
enddo
xyz_c = xyz_c / 4d0
- call lex_order(xyz_c,new_index_list,nnodes,f%nspec)
+ call lex_order(xyz_c,new_index_list,f%nspec)
f%ispec2 = f%ispec2(new_index_list)
f%inodes2 = f%inodes2(:,new_index_list)
end subroutine reorder_fault_elements_single
! ---------------------------------------------------------------------------------------------------
- subroutine lex_order(xyz_c,loc,nnodes,nspec)
+ subroutine lex_order(xyz_c,loc,nspec)
- integer, intent(in) :: nnodes,nspec
+ integer, intent(in) :: nspec
integer, intent(out) :: loc(nspec)
double precision, intent(in) :: xyz_c(3,nspec)
@@ -370,14 +368,14 @@
subroutine fault_repartition_not_parallel (nelmnts, nnodes, elmnts, nsize, nproc, part, esize)
- integer(long), intent(in) :: nelmnts,nsize
+ integer, intent(in) :: nelmnts,nsize
integer, intent(in) :: nnodes, nproc, esize
integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
integer, dimension(0:nelmnts-1), intent(inout) :: part
integer, dimension(0:nnodes-1) :: nnodes_elmnts
integer, dimension(0:nsize*nnodes-1) :: nodes_elmnts
- integer :: i,j, ipart,nproc_null,nproc_null_final
+ integer :: i,ipart,nproc_null,nproc_null_final
integer :: k1, k2, k,e,iflt,inode
integer, dimension(:), allocatable :: elem_proc_null
@@ -450,17 +448,16 @@
end subroutine fault_repartition_not_parallel
! ---------------------------------------------------------------------------------------------------
- subroutine fault_repartition_parallel (nelmnts, part, nodes_coords,nnodes)
+ subroutine fault_repartition_parallel (nelmnts, part, nodes_coords)
integer(long), intent(in) :: nelmnts
integer, dimension(0:nelmnts-1), intent(inout) :: part
- integer, intent(in) :: nnodes
double precision, dimension(3,nnodes), intent(in) :: nodes_coords
integer :: i,iflt,e,e1,e2,proc1,proc2
! Reorder both fault sides so that elements facing each other have the same index
- call reorder_fault_elements(nodes_coords,nnodes)
+ call reorder_fault_elements(nodes_coords)
!JPA loop over all faults
!JPA loop over all fault element pairs
@@ -495,7 +492,7 @@
integer, intent(in) :: IIN_database
integer, intent(in) :: iproc
- integer(long), intent(in) :: nelmnts
+ integer, intent(in) :: nelmnts
integer, dimension(0:nelmnts-1), intent(in) :: part
integer, dimension(0:nelmnts-1), intent(in) :: glob2loc_elmnts
integer, dimension(0:), intent(in) :: glob2loc_nodes_nparts
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_common.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_common.f90 2012-12-12 01:57:36 UTC (rev 21134)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_common.f90 2012-12-12 02:03:01 UTC (rev 21135)
@@ -42,14 +42,13 @@
!---------------------------------------------------------------------
-subroutine initialize_fault (bc,IIN_BIN,dt_tmp)
+subroutine initialize_fault (bc,IIN_BIN)
use specfem_par
use specfem_par_elastic, only : rmassx,rmassy,rmassz
class(fault_type), intent(inout) :: bc
integer, intent(in) :: IIN_BIN
- real(kind=CUSTOM_REAL), intent(in) :: dt_tmp
real(kind=CUSTOM_REAL) :: tmp_vec(3,NGLOB_AB)
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: jacobian2Dw
@@ -446,11 +445,11 @@
integer :: i,k,IOUT
character(len=10) :: my_fmt
- 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
+ integer, dimension(8) :: time_values
+ call date_and_time(VALUES=time_values)
+
IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
write(my_fmt,'(a,i1,a)') '(',dataT%ndat+1,'(E15.7))'
@@ -459,7 +458,7 @@
open(IOUT,file='../OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
write(IOUT,*) "# problem=TPV104" ! WARNING: this should be a user input
write(IOUT,*) "# author=Surendra Nadh Somala" ! WARNING: this should be a user input
- write(IOUT,1000) today(2), today(1), today(3), now
+ write(IOUT,1000) time_values(2), time_values(3), time_values(1), time_values(5), time_values(6), time_values(7)
write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
write(IOUT,*) "# code_version=1.1"
write(IOUT,*) "# element_size=100 m (*5 GLL nodes)" ! WARNING: this should be a user input
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_dynamic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_dynamic.f90 2012-12-12 01:57:36 UTC (rev 21134)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_dynamic.f90 2012-12-12 02:03:01 UTC (rev 21135)
@@ -198,7 +198,7 @@
NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
- call initialize_fault(bc,IIN_BIN,dt)
+ call initialize_fault(bc,IIN_BIN)
if (bc%nspec>0) then
@@ -427,7 +427,7 @@
real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T,dD,dV,dA
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength,tStick,tnew, &
theta_old, theta_new, dc, &
- ta,Vf_old,Vf_new,TxExt
+ Vf_old,Vf_new,TxExt
real(kind=CUSTOM_REAL) :: half_dt,TLoad,DTau0,GLoad,time
integer :: i
@@ -692,17 +692,16 @@
!=====================================================================
-subroutine rsf_init(f,T0,V,vel,nucFload,coord,IIN_PAR)
+subroutine rsf_init(f,T0,V,nucFload,coord,IIN_PAR)
type(rsf_type), intent(out) :: f
real(kind=CUSTOM_REAL), intent(in) :: T0(:,:)
real(kind=CUSTOM_REAL), intent(inout) :: V(:,:)
- real(kind=CUSTOM_REAL), intent(inout) :: vel(:,:)
real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
real(kind=CUSTOM_REAL), pointer :: nucFload(:)
integer, intent(in) :: IIN_PAR
- real(kind=CUSTOM_REAL) :: V0,f0,a,b,L,theta,theta_init,V_init,fw,Vw, C,T
+ real(kind=CUSTOM_REAL) :: V0,f0,a,b,L,theta_init,V_init,fw,Vw, C,T
integer :: nV0,nf0,na,nb,nL,nV_init,ntheta_init,nfw,nVw, nC,nForcedRup
real(kind=CUSTOM_REAL) :: W1,W2,w,hypo_z
real(kind=CUSTOM_REAL) :: x,z
@@ -872,7 +871,7 @@
type(rsf_type), intent(in) :: f
real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
real(kind=CUSTOM_REAL) :: mu(size(V))
- mu = f%a * asinh( V/2.0/f%V0 * exp((f%f0 + f%b*log(f%theta*f%V0/f%L))/f%a ) ) ! Regularized
+ mu = f%a * asinh_slatec( V/2.0/f%V0 * exp((f%f0 + f%b*log(f%theta*f%V0/f%L))/f%a ) ) ! Regularized
end function rsf_mu
@@ -888,7 +887,7 @@
else
arg = exp(theta/a)/TWO/V0
endif
- fn = tStick - Z*x - a*Seff*asinh(x*arg)
+ fn = tStick - Z*x - a*Seff*asinh_slatec(x*arg)
df = -Z - a*Seff/sqrt(ONE + (x*arg)**2)*arg
end subroutine funcd
@@ -991,19 +990,17 @@
!===============================================================
! OUTPUTS
-subroutine SCEC_Write_RuptureTime(dataXZ,DT,NT,iflt)
+subroutine SCEC_Write_RuptureTime(dataXZ,iflt)
type(dataXZ_type), intent(in) :: dataXZ
- real(kind=CUSTOM_REAL), intent(in) :: DT
- integer, intent(in) :: NT,iflt
+ integer, intent(in) :: iflt
integer :: i,IOUT
character(len=70) :: filename
- integer*4 today(3), now(3)
- call idate(today) ! today(1)=day, (2)=month, (3)=year
- call itime(now) ! now(1)=hour, (2)=minute, (3)=second
+ integer, dimension(8) :: time_values
+ call date_and_time(VALUES=time_values)
write(filename,"('../OUTPUT_FILES/RuptureTime_Fault',I0)") iflt
@@ -1012,7 +1009,7 @@
open(IOUT,file=trim(filename),status='replace')
write(IOUT,*) "# problem=TPV104"
write(IOUT,*) "# author=Surendra Nadh Somala"
- write(IOUT,1000) today(2), today(1), today(3), now
+ write(IOUT,1000) time_values(2), time_values(3), time_values(1), time_values(5), time_values(6), time_values(7)
write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
write(IOUT,*) "# code_version=1.1"
write(IOUT,*) "# element_size=100 m (*5 GLL nodes)"
@@ -1194,5 +1191,159 @@
end subroutine write_dataXZ
+!---------------------------------------------------------------
+
+! asinh() function taken from Netlib
+! April 1977 edition. W. Fullerton, C3, Los Alamos Scientific Lab.
+
+! taken from http://www.tddft.org/trac/octopus/browser/trunk/src/asinh.F90?rev=2
+
+! and modified by Dimitri Komatitsch in December 2012 for portability
+
+ double precision function asinh_slatec(x)
+
+ double precision, intent(in) :: x
+
+ integer, parameter :: NSERIES = 39
+
+ double precision, parameter :: asnhcs(NSERIES) = (/ &
+ -.12820039911738186343372127359268D+0, -.58811761189951767565211757138362D-1, &
+ +.47274654322124815640725249756029D-2, -.49383631626536172101360174790273D-3, &
+ +.58506207058557412287494835259321D-4, -.74669983289313681354755069217188D-5, &
+ +.10011693583558199265966192015812D-5, -.13903543858708333608616472258886D-6, &
+ +.19823169483172793547317360237148D-7, -.28847468417848843612747272800317D-8, &
+ +.42672965467159937953457514995907D-9, -.63976084654366357868752632309681D-10, &
+ +.96991686089064704147878293131179D-11, -.14844276972043770830246658365696D-11, &
+ +.22903737939027447988040184378983D-12, -.35588395132732645159978942651310D-13, &
+ +.55639694080056789953374539088554D-14, -.87462509599624678045666593520162D-15, &
+ +.13815248844526692155868802298129D-15, -.21916688282900363984955142264149D-16, &
+ +.34904658524827565638313923706880D-17, -.55785788400895742439630157032106D-18, &
+ +.89445146617134012551050882798933D-19, -.14383426346571317305551845239466D-19, &
+ +.23191811872169963036326144682666D-20, -.37487007953314343674570604543999D-21, &
+ +.60732109822064279404549242880000D-22, -.98599402764633583177370173440000D-23, &
+ +.16039217452788496315232638293333D-23, -.26138847350287686596716134399999D-24, &
+ +.42670849606857390833358165333333D-25, -.69770217039185243299730773333333D-26, &
+ +.11425088336806858659812693333333D-26, -.18735292078860968933021013333333D-27, &
+ +.30763584414464922794065920000000D-28, -.50577364031639824787046399999999D-29, &
+ +.83250754712689142224213333333333D-30, -.13718457282501044163925333333333D-30, &
+ +.22629868426552784104106666666666D-31 /)
+
+ double precision, parameter :: aln2 = 0.69314718055994530941723212145818D0
+
+! series for asnh on the interval 0. to 1.00000d+00
+! with weighted error 2.19e-17
+! log weighted error 16.66
+! significant figures required 15.60
+! decimal places required 17.31
+!
+
+ integer, save :: nterms = 0
+ double precision, save :: xmax = 0.d0, sqeps = 0.d0
+
+! taken from http://people.sc.fsu.edu/~jburkardt/f_src/machine/machine.f90
+ double precision, parameter :: d1mach_3 = 1.110223024625157D-016
+
+ integer, external :: inits
+ double precision, external :: csevl
+ double precision :: y
+
+ if (nterms == 0) then
+ nterms = inits(asnhcs, NSERIES, 0.1d0*d1mach_3)
+ sqeps = sqrt(d1mach_3)
+ xmax = 1.d0/sqeps
+ endif
+
+ y = abs(x)
+ if (y <= 1.d0) then
+ asinh_slatec = x
+ if (y > sqeps) asinh_slatec = x*(1.d0 + csevl(2.d0*x*x-1.d0, asnhcs, nterms))
+ return
+ endif
+
+ if (y < xmax) asinh_slatec = log(y + sqrt(y**2 + 1.d0))
+ if (y >= xmax) asinh_slatec = aln2 + log(y)
+ asinh_slatec = sign(asinh_slatec, x)
+
+end function asinh_slatec
+
+
+! April 1977 version. W. Fullerton, C3, Los Alamos Scientific Lab.
+! Evaluate the n-term Chebyshev series cs at x. Adapted from
+! R. Broucke, Algorithm 446, C.A.C.M., 16, 254 (1973). Also see Fox
+! and Parker, Chebyshev polynomials in numerical analysis, Oxford Press, p.56.
+!
+! input arguments --
+! x value at which the series is to be evaluated.
+! cs array of n terms of a Chebyshev series.
+! in evaluating cs, only half the first coefficient is summed.
+! n number of terms in array cs.
+
+double precision function csevl(x, cs, n)
+
+ integer, intent(in) :: n
+ double precision, intent(in) :: x
+ double precision, intent(in) :: cs(n)
+
+ integer i, ni
+ double precision :: b0, b1, b2, twox
+
+ if (n < 1) stop 'Math::csevl: number of terms <= 0'
+ if (n > 1000) stop 'Math::csevl: number of terms > 1000'
+
+ if (x < -1.1d0 .or. x > 1.1d0) stop 'Math::csevl: x outside (-1,+1)'
+
+ b1 = 0.d0
+ b0 = 0.d0
+ twox = 2.d0*x
+
+ do i = 1, n
+ b2 = b1
+ b1 = b0
+ ni = n + 1 - i
+ b0 = twox*b1 - b2 + cs(ni)
+ enddo
+
+ csevl = 0.5d0 * (b0 - b2)
+
+end function csevl
+
+
+! April 1977 version. W. Fullerton, C3, Los Alamos Scientific Lab.
+!
+! Initialize the orthogonal series so that inits is the number of terms
+! needed to ensure that the error is no larger than eta. Ordinarily, eta
+! will be chosen to be one-tenth machine precision.
+!
+! input arguments --
+! os array of nos coefficients in an orthogonal series.
+! nos number of coefficients in os.
+! eta requested accuracy of series.
+
+integer function inits(os, nos, eta)
+
+ integer, intent(in) :: nos
+ double precision, intent(in) :: os(nos)
+ double precision, intent(in) :: eta
+
+ integer :: i, ii
+ double precision :: err
+
+ if (nos < 1) stop 'Math::inits: number of terms <= 0'
+
+ err = 0.d0
+ do ii=1,nos
+ i = nos + 1 - ii
+ err = err + abs(os(i))
+ if (err > eta) exit
+ enddo
+
+!!!!!!! if (i == nos) print *,'warning: Math::inits: eta may be too small'
+
+ inits = i
+
+end function inits
+
+
+
end module fault_solver_dynamic
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_kinematic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_kinematic.f90 2012-12-12 01:57:36 UTC (rev 21134)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_kinematic.f90 2012-12-12 02:03:01 UTC (rev 21135)
@@ -138,7 +138,7 @@
NAMELIST / KINPAR / kindt
- call initialize_fault(bc,IIN_BIN,dt)
+ call initialize_fault(bc,IIN_BIN)
if (bc%nspec>0) then
@@ -335,9 +335,9 @@
! OUTPUT v : slip rate on receivers.
-subroutine load_vslip_snapshots(dataXZ,itime,nglob,iflt)
+subroutine load_vslip_snapshots(dataXZ,itime,iflt)
- integer, intent(in) :: itime,nglob,iflt
+ integer, intent(in) :: itime,iflt
type(dataXZ_type), intent(inout) :: dataXZ
character(len=70) :: filename
integer :: IIN_BIN,ier,IOUT
More information about the CIG-COMMITS
mailing list