王培杰 发表于 2014-4-20 19:53:53

MT一维正演,用SOR法不收敛

原理如图,结果一直不收敛
代码如下program main
implicit none
integer,parameter::n=1000                     !地上n层 地下n层
complex*16,parameter::i=(0d0,1d0)               
real*16,parameter::PI=3.1415926535897932384626
real*16,parameter::miu=4d0*PI*1E-7            
real*16,parameter::w=1.5                        !SOR迭代因子
integer::j,k
real*16::o(-n-1:n),o1(-n-1:n),dz(-n-1:n+1),s,f,E14,ro,fa
!o地层电导率,o1分界面电导率,dz层厚度,f频率,E14四分之一厚度处电场,ro视电阻率,fa相位
complex*16::E(-n-1:n+1),H(-n-1:n),E1(-n-1:n+1),Hs,a,b
!E电场,H磁场,Hs地表磁场
open(123,file='a.dat')
dz=1d0                  !单层厚度1m
dz(-n-1)=1.0E6            !最上层厚度1000km,视作无限远
dz(n)=1.0E6               !最底层厚度1000km,视作无限远
forall (k=-n-1:-1:1) o(k)=1E-10      !空气中电导率1E-10,视作无限大
forall(k=0:300:1) o(k)=0.1d0         !0到300米,电导率0.1
forall(k=301:600:1) o(k)=0.005d0   !301到600米,电导率0.005
forall(k=601:n:1) o(k)=0.1d0         !601到无穷远,电导率0.1
do k=-n,n
    o1(k)=(o(k)*dz(k)+o(k-1)*dz(k-1))/(dz(k)+dz(k-1))!计算层界面电导率
enddo
do j=20,2000,20                     !j为厚度
    f=(500**2)*100/j                  !趋附深度求取频率
    E=(1d0,1d0)                     !迭代初值
    E(-n-1)=(1.0E-6,0d0)            !地上无穷远处电场1E-6v/m
    E(n+1)=(0d0,0d0)                !地下无穷远处电场为0
    do                              !SOR迭代法
      E1=E
      s=0
      do k=-n,n
            a=2d0*(E(k+1)/(dz(k)*dz(k)+dz(k-1)*dz(k))+E(k-1)/(dz(k-1)*dz(k)+dz(k-1)*dz(k-1)))
            b=2d0/(dz(k)*dz(k)+dz(k-1)*dz(k))+2d0/(dz(k-1)*dz(k)+dz(k-1)*dz(k-1))+i*f*miu*o1(k)
            E(k)=(1-w)*E(k)+w*a/b   
            s=s+abs(E(k)-E1(k))
      enddo
      if(s<1E-20)exit             !迭代结束条件
    enddo
    do k=-n-1,n
      H(k)=(E(k+1)-E(k))/(dz(k)*i*f*miu)      !计算各地层中点磁场
    enddo
    E14=0.75*E(0)+0.25d0*E(1)
    Hs=E14*o1(0)*dz(0)/2d0+H(0)               !地表处磁场
    ro=(abs(E(0)/Hs))**2/f/miu                  !卡尼亚视电阻率
    write(*,*)f,ro
    write(123,*)f,ro
enddo
pause
endprogram

aliouying 发表于 2014-4-21 10:21:50

本帖最后由 aliouying 于 2014-4-21 10:23 编辑

为什么要用差分呢? sor迭代收敛很慢,一般需要迭代1000次以上,甚至几千次,你的程序有些问题,大致看了下,改了改program main
implicit none
integer,parameter::n=1000                     !地上n层 地下n层
complex*16,parameter::i=(0d0,1d0)               
real*8,parameter::PI=3.1415926535897932384626d0
real*8,parameter::miu=4.0d0*PI*1.0d-7            
real*8,parameter::w=0.8                        !SOR迭代因子
integer::j,k
real*8::o(-n-1:n),o1(-n:n),dz(-n-1:n+1),s,f,E14,ro,fa
!o地层电导率,o1分界面电导率,dz层厚度,f频率,E14四分之一厚度处电场,ro视电阻率,fa相位
complex*16::E(-n-1:n+1),H(-n-1:n),E1(-n-1:n+1),Hs,a,b
!E电场,H磁场,Hs地表磁场
open(123,file='a.dat')
dz=1d0                  !单层厚度1m
dz(-n-1)=1.0E6            !最上层厚度1000km,视作无限远
dz(n)=1.0E6               !最底层厚度1000km,视作无限远
forall (k=-n-1:-1:1) o(k)=1E-10      !空气中电导率1E-10,视作无限大
forall(k=0:300:1) o(k)=0.1d0         !0到300米,电导率0.1
forall(k=301:600:1) o(k)=0.005d0   !301到600米,电导率0.005
forall(k=601:n:1) o(k)=0.1d0         !601到无穷远,电导率0.1
do k=-n,n
    o1(k)=(o(k)*dz(k)+o(k-1)*dz(k-1))/(dz(k)+dz(k-1))!计算层界面电导率
enddo
do j=20,2000,20                     !j为厚度
    f=2*PI*(500**2)*100/j         !趋附深度求取频率
    E=(1d0,1d0)                     !迭代初值
    !E(-n-1)=(1.0E-6,0d0)         !地上无穷远处电场1E-6v/m
    E(n+1)=(0d0,0d0)                !地下无穷远处电场为0
    do                              !SOR迭代法
      E1=E
      s=0
      do k=-n,n
            a=2d0*(E(k+1)/(dz(k)*dz(k)+dz(k-1)*dz(k))+E(k-1)/(dz(k-1)*dz(k)+dz(k-1)*dz(k-1)))
            b=2d0/(dz(k)*dz(k)+dz(k-1)*dz(k))+2d0/(dz(k-1)*dz(k)+dz(k-1)*dz(k-1))+i*f*miu*o1(k)
            E(k)=(1-w)*E(k)+w*a/b   
            s=s+abs(E(k)-E1(k))**2
      enddo
      s = dsqrt( s/n/2 )
      if(s<1E-10)exit             !迭代结束条件
    enddo
    do k=-n-1,n
      H(k)=(E(k+1)-E(k))/(dz(k)*i*f*miu)      !计算各地层中点磁场
    enddo
    E14=0.75*E(0)+0.25d0*E(1)
    Hs=E14*o1(0)*dz(0)/2d0+H(0)               !地表处磁场
    ro=(abs(E(0)/Hs))**2/f/miu                  !卡尼亚视电阻率
    write(*,*)f,ro
    write(123,*)f,ro
enddo
pause
endprogram



geophylika 发表于 2014-5-20 00:04:00

做物探的不少,其实f 的代码现成的不少`
页: [1]
查看完整版本: MT一维正演,用SOR法不收敛