Fortran Coder

标题: 三种有限元差分法求散度 [打印本页]
作者: 云卷云舒    时间: 2017-12-6 20:11
标题: 三种有限元差分法求散度
本帖最后由 云卷云舒 于 2017-12-6 20:11 编辑

各位大神，新手求助求解散度！有经向纬向两组速度数据，（附加图已给出,这是我读出数据后又整理了一下输出的，在我读取源文件是是先读的一行x）
用前差、后差、中心差求解各点的散度（偏u偏x+偏v偏y）（由于边界没法求，有的点没有散度）
但是我求出来三种方法的散度值完全不一样，不知道问题出在那里，求高手指点一下，万分感激！！！
一下是我的程序：（只看计算部分即可，print部分可以不看）

program test11
   implicit none
   integer::n,i,j
   real::x(17),y(17)
   real,parameter::d=0.25
   real::u(17,17),v(17,17)
   real::df_ux(16,16),df_vy(16,16),df(16,16)
   real::da_ux(2:17,2:17),da_vy(2:17,2:17),da(2:17,2:17)
   real::dm_ux(2:16,2:16),dm_vy(2:16,2:16),dm(2:16,2:16)!----------------一种方法不能只开一个数组，因为两次求导的内外层循环顺序反了！！！！
   open(unit=10,file='u.txt',status='old')
   open(unit=20,file='v.txt',status='old')
   read(10,*) (x(i),i=1,17)  !-----------注意这里
   do i=1,17
       read(20,'(18f12.7)') y(i),v(i,:)
       read(10,'(18f12.7)') y(i),u(i,:)      
   end do
   !============================计算====================================
   !=============================前差=======================
   !---------u/x
    do j=1,16
       do i=1,16
           df_ux(i,j)=(u(i+1,j)-u(i,j))/d!----i对应x,j对应y
       end do
   end do
   !---------v/y
   do i=1,16
       do j=1,16
           df_vy(i,j)=(v(i,j+1)-v(i,j))/d
       end do
   end do
     df=df_ux+df_vy
   !=============================后差=======================
    !---------u/x
   do j=2,17
       do i=17,2,-1
           da_ux(i,j)=(u(i-1,j)-u(i,j))/d
       end do
   end do
   !---------v/y
   do i=2,17
       do j=17,2,-1
           da_vy(i,j)=(v(i,j-1)-v(i,j))/d
       end do
   end do
     da=da_ux+da_vy
   !===========================中心差=======================
   !---------u/x
   do j=2,16
       do i=2,16
           dm_ux(i,j)=(u(i+1,j)-u(i-1,j))/(2.0*d)
       end do
   end do
   !---------v/y
   do j=2,16
       do i=2,16
           dm_vy(i,j)=(v(i,j+1)-v(i,j-1))/(2.0*d)
       end do
   end do
     dm=dm_ux+dm_vy

   !==========================print===========================================
   !=============================前差=======================
   open(unit=11,file='div_f.txt')

   do i=16,1,-1
       write(11,'(18f12.7)') y(i),df(i,:)  
   end do
   write(11,'(12x,16f12.7)') (x(i),i=1,16)
   !-------------------------------------
   print*,"前差："

   do i=16,1,-1
       print'(18f10.7)', y(i),df(i,:)   
   end do
   print'(10x,16f10.7)',(x(i),i=1,16)
   !=============================后差=======================
   open(unit=12,file='div_a.txt')
   do i=17,2,-1
       write(12,'(18f12.7)')y(i),da(i,:)   
   end do
   write(12,'(12x,16f12.7)')(x(i),i=2,17)
   !---------------------------------------
     print*,"后差："

   do i=17,2,-1
       print'(18f10.7)', y(i),da(i,:)   
   end do

   print'(10x,16f10.7)',(x(i),i=2,17)
   !===========================中心差=======================
   open(unit=13,file='div_m.txt')
   do i=16,2,-1
       write(13,'(18f12.7)') y(i),dm(i,:)  
   end do
   write(13,'(12x,15f12.7)') (x(i),i=2,16)
    !-------------------------------
     print*,"中心差："

   do i=16,2,-1
       print'(18f10.7)', y(i),dm(i,:)   
   end do
   print'(10x,15f10.7)', (x(i),i=2,16)
   end program


97f899cffb375c2aaf3102c4af5c791.png (53.55 KB, 下载次数: 264)

以u为例，左侧是y,最后一行是x

以u为例，左侧是y,最后一行是x
作者: li913    时间: 2017-12-6 21:44
各种差分格式的结果有差异，这是肯定的，不然就没必要那么多格式了。
后差写错了，符号不对。



欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2