zhangyiss/DSurfTomo
1
0
Fork 0
mirror of https://github.com/HongjianFang/DSurfTomo.git synced 2025-05-07 15:41:14 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

clean version of main.f90

Browse Source
This commit is contained in:
Hongjian Fang 2019-01-15 08:57:51 -05:00
parent 404f55cd0f
commit 82ebf4283d
2 changed files with 556 additions and 713 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
src/CalSurfG.f90
View File

@ -940,7 +940,7 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
goxdf,gozdf,dvxdf,dvzdf,kmaxRc,kmaxRg,kmaxLc,kmaxLg, &
tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk, &
scxf,sczf,rcxf,rczf,nrc1,nsrcsurf1,kmax,nsrcsurf,nrcf, &
nar,writepath)
nar)
USE globalp
USE traveltime
IMPLICIT NONE
@ -1027,7 +1027,6 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
integer ii,jj,kk,nn,istep
integer level,maxlevel,maxleveld,HorizonType,VerticalType,PorS
real,parameter::ftol=1e-4
integer writepath
integer ig, igroup
gdx=8
@ -1382,7 +1381,7 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
call gridder(velf0)
endif
count11=count11+1
CALL rpaths(x,z,fdm,rcxf(istep,srcnum,knumi),rczf(istep,srcnum,knumi),writepath)
CALL rpaths(x,z,fdm,rcxf(istep,srcnum,knumi),rczf(istep,srcnum,knumi))
row(1:nparpi)=0.0
do jj=1,nvz
do kk=1,nvx
@ -1747,7 +1746,7 @@ END SUBROUTINE srtimes
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!SUBROUTINE rpaths(wrgf,csid,cfd,scx,scz)
!SUBROUTINE rpaths()
SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz,writepath)
SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz)
USE globalp
IMPLICIT NONE
INTEGER, PARAMETER :: i5=SELECTED_REAL_KIND(5,10)
@ -1768,7 +1767,6 @@ SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz,writepath)
!fang!------------------------------------------------
real fdm(0:nvz+1,0:nvx+1)
REAL(KIND=i10) surfrcx,surfrcz
integer writepath
!fang!------------------------------------------------
!
! ipx,ipz = Coordinates of cell containing current point
@ -2253,14 +2251,14 @@ SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz,writepath)
! Write ray paths to output file
!
!fang! IF(wrgf.EQ.csid.OR.wrgf.LT.0)THEN
if(writepath == 1) then
WRITE(40,*)'#',nrp
DO j=1,nrp
rayx=(pi/2-rgx(j))*180.0/pi
rayz=rgz(j)*180.0/pi
WRITE(40,*)rayx,rayz
ENDDO
endif
!if(writepath == 1) then
! WRITE(40,*)'#',nrp
! DO j=1,nrp
! rayx=(pi/2-rgx(j))*180.0/pi
! rayz=rgz(j)*180.0/pi
! WRITE(40,*)rayx,rayz
! ENDDO
!endif
!fang! ENDIF
!
! Write partial derivatives to output file

157
src/main.f90
View File

@ -96,26 +96,12 @@ program SurfTomo
real spfra
real noiselevel
integer ifsyn
integer writepath
real averdws
real maxnorm
real threshold0
! FOR MODEL VARIATION
!------------------------------------------------
integer idx
integer counte
real stdvs
! integer numrand
! real,allocatable,dimension(:,:)::modstat
! real,allocatable,dimension(:)::modsig
real gaussian
external gaussian
integer modest
counte=0
! OPEN FILES FIRST TO OUTPUT THE PROCESS
open(34,file='IterVel.out')
nout=36
open(nout,file='lsmr.txt')
@ -218,8 +204,6 @@ program SurfTomo
read(10,*)ifsyn
read(10,*)noiselevel
read(10,*) threshold0
! read(10,*) modest
! read(10,*) numrand
close(10)
nrc=nsrc
kmax=kmaxRc+kmaxRg+kmaxLc+kmaxLg
@ -301,10 +285,6 @@ program SurfTomo
allocate(norm(maxvp), stat=checkstat)
allocate(vsf(nx,ny,nz), stat=checkstat)
allocate(vsftrue(nx,ny,nz), stat=checkstat)
! FOR MODEL VARIATION
!------------------------------------------------
! allocate(modstat(numrand,maxvp))
! allocate(modsig(maxvp))
allocate(rw(maxnar), stat=checkstat)
if(checkstat > 0)then
@ -335,8 +315,6 @@ program SurfTomo
write(*,*) 'grid points in depth direction:(km)'
write(*,'(50f7.1)') depz
! CHECKERBOARD TEST
if (ifsyn == 1) then
write(*,*) 'Checkerboard Resolution Test Begin'
@ -360,57 +338,24 @@ program SurfTomo
! ITERATE UNTILL CONVERGE
writepath = 0
do iter = 1,maxiter
iw = 0
rw = 0.0
col = 0
! COMPUTE SENSITIVITY MATRIX
if (iter == maxiter) then
writepath = 1
open(40,file='raypath.out')
endif
write(*,*) 'computing sensitivity matrix...'
call CalSurfG(nx,ny,nz,maxvp,vsf,iw,rw,col,dsyn,&
goxd,gozd,dvxd,dvzd,kmaxRc,kmaxRg,kmaxLc,kmaxLg,&
tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk,&
scxf,sczf,rcxf,rczf,nrc1,nsrc1,kmax,&
nsrc,nrc,nar,writepath)
nsrc,nrc,nar)
do i = 1,dall
cbst(i) = obst(i) - dsyn(i)
enddo
! write out rw, iw and cbst for testing
!open(44,file='iw.dat')
!write(44,*) nar
!do i = 2,nar+1
!write(44,*) iw(i)
!enddo
!do i = 1,nar
!write(44,*) col(i)
!enddo
!close(44)
!open(44,file='rw.dat')
!do i = 1,nar
!write(44,*) rw(i)
!enddo
!close(44)
!open(44,file='residal.dat')
!do i = 1,dall
!write(44,*) cbst(i)
!enddo
!close(44)
!threshold = threshold0+(maxiter/2-iter)/3*0.5
do i = 1,dall
! datweight(i) = 1.0
! if(abs(cbst(i)) > threshold) then
! datweight(i) = exp(-(abs(cbst(i))-threshold))
! endif
datweight(i) = 0.01+1.0/(1+0.05*exp(cbst(i)**2*threshold0))
cbst(i) = cbst(i)*datweight(i)
enddo
@ -431,7 +376,6 @@ program SurfTomo
enddo
averdws=averdws/maxvp
write(66,*)'Maximum and Average DWS values:',maxnorm,averdws
!write(66,*)'Threshold is:',threshold
! WRITE OUT RESIDUAL FOR THE FIRST AND LAST ITERATION
if(iter.eq.1) then
@ -453,7 +397,6 @@ program SurfTomo
! ADDING REGULARIZATION TERM
!weight=dnrm2(dall,cbst,1)**2/dall*weight0
weight=weight0
nar_tmp=nar
nars=0
@ -528,7 +471,6 @@ program SurfTomo
call LSMR(m, n, leniw, lenrw,iw,rw,cbst, damp,&
atol, btol, conlim, itnlim, localSize, nout,&
dv, istop, itn, anorm, acond, rnorm, arnorm, xnorm)
!if(istop==3) print*,'istop = 3, large condition number'
mean = sum(cbst(1:dall))/dall
std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2)
@ -544,8 +486,6 @@ program SurfTomo
mean = sum(cbst(1:dall))/dall
std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2)
!write(*,'(i2,a)'),iter,'th iteration...'
!write(*,'(a,f7.3)'),'weight is:',weight
write(*,'(a,f8.1,a,f8.2,a,f8.3)'),'residual before weighting: ',mean*1000,'ms ',1000*std_devs,'ms ',&
dnrm2(dall,cbst,1)/sqrt(real(dall))
write(66,'(i2,a)'),iter,'th iteration...'
@ -574,18 +514,6 @@ program SurfTomo
enddo
enddo
enddo
write(34,*)'OUTPUT S VELOCITY AT ITERATION',iter
do k=1,nz
do j=1,ny
write(34,'(100f7.3)') (vsf(i,j,k),i=1,nx)
enddo
enddo
write(34,*)',OUTPUT DWS AT ITERATION',iter
do k=1,nz-1
do j=2,ny-1
write(34,'(100f8.1)') (norm((k-1)*(ny-2)*(nx-2)+(j-2)*(nx-2)+i-1),i=2,nx-1)
enddo
enddo
write(outmodel,'(a,a,i3.3)') trim(inputfile),'Measure.dat.iter',iter
open(64,file=outmodel)
@ -645,77 +573,10 @@ program SurfTomo
to ',outmodel
endif
close(34)
close(40)
close(nout) !close lsmr.txt
close(66) !close surf_tomo.log
!! USE RANDOM MODEL TO OBTAIN THE MODEL VARIATION
! !modest = 1
! if (modest ==1) then
!
! write(*,*) 'model variation estimation begin...'
! do iter = 1,numrand
! call init_random_seed()
! vsftrue=vsf
! DO K=1,NZ-1
! DO J=2,NY-1
! DO I=2,NX-1
! idx = (k-1)*(ny-2)*(nx-2)+(j-2)*(nx-2)+i-1
! dv(idx) = 0.1/EXP(2*NORM(idx)/maxnorm)*gaussian()
! VSFTRUE(I,J,K) = VSF(I,J,K)+dv(idx)
! ENDDO
! ENDDO
! ENDDO
! write(*,*),'maximum and minimum velocity variation',maxval(dv),minval(dv)
!
! call synthetic(nx,ny,nz,maxvp,vsftrue,dsyn,&
! goxd,gozd,dvxd,dvzd,kmaxRc,kmaxRg,kmaxLc,kmaxLg,&
! tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk,&
! scxf,sczf,rcxf,rczf,nrc1,nsrc1,kmax,&
! nsrc,nrc,0.0)
!
! do i = 1,dall
! cbst(i) = obst(i) - dsyn(i)
! enddo
!
! write(*,*), dnrm2(dall,cbst,1)/sqrt(real(dall)), 1.05*std_devs
! if (dnrm2(dall,cbst,1)/sqrt(real(dall)) < 1.05*std_devs) then
! counte = counte + 1
! modstat(counte,:) = dv
! endif
!
! enddo ! iteration for random models
!
! write(*,*),'number of of models satisfy requirements',counte
! modsig = 1.0
! if (counte>0) then
! do i=1,maxvp
! !statis
! !mean = sum(cbst(1:dall))/dall
! !std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2)
! mean = sum(modstat(1:counte,i))/counte
! stdvs = sqrt(sum(modstat(1:counte,i)**2)/counte-mean**2)
! modsig(i) = stdvs
! enddo
! endif
!
! write(*,*),'write model variation to "model_variation.dat"'
! open (64,file='model_variation.dat')
! do k=1,nz-1
! do j=1,ny-2
! do i=1,nx-2
! idx = (k-1)*(ny-2)*(nx-2)+(j-1)*(nx-2)+i
! write(64,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),modsig(idx)
! enddo
! enddo
! enddo
! close(64)
! write(*,*) 'finishing model variation estimation'
! endif
deallocate(obst)
deallocate(dsyn)
deallocate(dist)
@ -746,19 +607,3 @@ program SurfTomo
end program
!-----------------------------------------------------------------------
! Generate seed for random number generator of fortran
! Note: only need to be called once inside one program
!-----------------------------------------------------------------------
subroutine init_random_seed()
integer :: i,n,clock
integer,dimension(:),allocatable :: seed
call random_seed(size=n)
allocate(seed(n))
call system_clock(count=clock)
seed=clock+37*(/(i-1,i=1,n)/)
call random_seed(PUT=seed)
deallocate(seed)
end subroutine