作者Gill0920 (不斷邁進!!)
看板Fortran
標題[問題]麻煩幫我看看~~非常感謝!!!
時間Sun Feb 28 22:07:52 2010
這是我的程式 前面算係數的部分都確認過了
最後的三個大廻圈 不知道要怎麼表示跟輸出
以下是我的問題
1. Q1 一直是近似值，應該要是0.045 (但一直都是0.04499999 除了這個係數 其他的都對
而且不會有這種進位的問題 )
2. 把時間迴圈另外寫的表示方法? 我寫 BS(NT)--->不行!@@" 如何寫才正確?
3. 由於BS.B1.B2 皆為EFDM 由時刻t.t+1 求t+2 所以需要每秒輸出一次結果再進行取代
請問這樣的寫法順序上是否正確? (且輸出的位置正確嗎?)
4. error 測試是跟前一項的誤差在1%內就可輸出，如果沒有就重跑時間迴圈，但
執行出來的錯誤是說重複時間迴圈的動作了@@" 這樣的想法不對?!還是要連
空間上的迴圈一起重跑呢?
! 最後謝謝熱心幫忙的大大們~~ 祝大家都除蟲順利!!
PROGRAM EX
USE IMSL
REAL*8 A(3,3),B(3,3),X(3,3)
REAL*8 PS(401,301),P1(401,301),P2(401,301),BS(401,301),B1(401,301),B2(401,301)
!!
integer*8 kb,ks,g,k1,k2,T
!!
REAL*8 M,n1,n2,n3,kr1,kr2,pi,dt,dx,dy,f,y1
!!
PARAMETER( KI=401,JM=301,NM=100001,pi=3.14159 )
OPEN (1,FILE='SOLID.TXT',status='unknown')
OPEN (2,FILE='FLUID1.TXT',status='unknown')
OPEN (3,FILE='FLUID2.TXT',status='unknown')
OPEN (4,FILE='all.TXT',status='unknown')
WRITE(*,*)KI,JM,NM
DT =1D-5
DX=1D-2
DY=1D-2
OUTDT=1
write(*,*)DT,DX,DY
!!fitting coefficients E,M,Y,F,X1!!
E = 2.145
M = 1.-(1./E)
Y = 0.5
F= 53D-14
X1 = 1.0
write(*,*)E,M,F
!!elastic coefficients KB,KS,K1,K2,G!!
KB = 8.33E6
KS = 35.E9
K1 = 0.145E6
K2 = 2.25E9
G = 3.85E6
write(*,*)kb,ks,k1,k2,g
!!density PI1,PI2,PIS!!
PI1 = 1.1
PIS = 2650.0
PI2 = 997.0
write(*,*)PI1,PIS,PI2
!!viscosity Y1,Y2!!
Y1 = 18D-6
Y2 = 0.001
write(*,*)Y1,Y2
!!porosity H!!
H = 0.45
S1 =0.1
S2 =1.-S1
IF (S1==0) THEN
DPS=0
ELSE IF (S1==1.0) THEN
DPS=-1000000000000;
ELSE
DPS=(PI2*9.81/(M*E*X1))*(((1.-S1)**(-E/(E-1.))-1.)**((1.-E)/E))*((1.-S1)**(-(2.
*E-1)/(E-1.)));
END IF
write(*,*)H,DPS,S1,S2
!!compare with Tuncay & Corapcioglu!!
N1 = (KS*(1.-h))-KB
N2 = DPS*S1*(1.-S1)
N3 = (N1*((K1*N2*S1)+(K1*K2)+(K2*N2*(1.-S1))))+(1225.E18)*H*(K1*(1.-S1)+N2+(K2*
S1))
write(*,*)N1,N2,N3
!!symmetric matrix!!
U11 = KS/N3*((1.-H)*N1*(K1*K2+K1*N2*S1+K2*N2*(1.-S1))+KB*KS*H*(K1*(1.-S1)+K2*S1
+N2))
U12 = N1*KS*H*K1*S1*(K2+N2)/N3
U21 = N1*KS*H*K1*S1*(K2+N2)/N3
U13 = N1*KS*H*K2*(1.-S1)*(K1+N2)/N3
U31 = N1*KS*H*K2*(1.-S1)*(K1+N2)/N3
U22 = (K1*(S1**2)*H*(1225.E18)*H*K2+K1*(S1**2)*(H**2)*N2*(1225.E18)+K1*S1*H*K2*
N1*N2*(1.-S1))/N3
U23 = (-K1)*K2*H*S1*(1.0-S1)*(N1*N2-H*(1225.E18))/N3
U32 = (-K1)*K2*H*S1*(1.-S1)*(N1*N2-H*(1225.E18))/N3
U33 = (K1*K2*(1225.E18)*((1.-S1)**2)*(H**2)+(1225.E18)*N2*(H**2)*K2*(1.-S1)+K1*
S1*N1*N2*K2*(1.-S1)*H)/N3
write(*,*)U11,U12,U21,U13,U31,U22,U23,U32,U33
!!inertiaL coupling RS,AA11,AA12,AA22,AA21,A111!!
RS = 0.5*(1.+1./H)
AA11 = PI1*(S1)*H-RS*PI1*(S1)*H
AA12 = -0.1*SQRT(((RS)**2)*PI1*PI2*(S1*H)*(S2*H))
AA21 = -0.1*SQRT(((RS)**2)*PI1*PI2*(S1*H)*(S2*H))
AA22 = PI2*(1.-S1)*H-RS*PI2*(1.-S1)*H
write(*,*)RS,AA11,AA12,AA21,AA22
!!permeability KR1,KR2!!
KR1 = ((1.-S2)**Y)*((1.-(S2)**(1./M))**(2*M))
KR2 = ((S2)**Y)*((1.-((1.-(S2)**(1./M))**M))**2)
T1 = KR1/Y1
T2 = KR2/Y2
write(*,*)KR1,KR2,T1,T2
!!volume fractIon Q1,Q2,QS!!
Q1 = S1*H
Q2 = S2*H
QS = 1.-Q1-Q2
write(*,*)Q1,Q2,QS
!!drag coefficients R11,R22!!
IF (KR1 ==0) THEN
R11 = 0
ELSE
R11 = -((Q1)**2)/(F*T1)
END IF
IF (KR2 == 0) THEN
R22 = 0;
ELSE
R22 = -((Q2)**2)/(F*T2)
END IF
write(*,*)R11,R22
!!the elements from stress-strain relations!!
D =3*(U11+U12+U13)*U22*U33+3*(U21+U22+U23)*U13*U23+3*(U31+U32+U33)*U12*U23-3*
(U11+U12+U13)*(U23**2)-3*(U12+U22+U23)*U12*U33-3*(U13+U23+U33)*U22*U13
C11 = (U22*U33-U23**2)/D
C12 = (3*(U13+U23+U33)*U23-3*(U12+U22+U23)*U33)/D
C13 = (3*(U12+U22+U23)*U23-3*(U13+U23+U33)*U22)/D
C21 = (U13*U23-U12*U33)/D
C22 = (3*(U11+U12+U13)*U33-3*(U13+U23+U33)*U13)/D
C23 = (3*(U13+U23+U33)*U12-3*(U11+U12+U13)*U23)/D
C31 = (U12*U23-U13*U22)/D
C32 = (3*(U12+U22+U23)*U13-3*(U11+U12+U13)*U23)/D
C33 = (3*(U11+U12+U13)*U22-3*(U12+U22+U23)*U12)/D
write(*,*)D,C11,C12,C13,C21,C22,C23,C31,C32,C33
!!T+1 term coefficients:change !!
B11 =(((PIS*QS)-AA11-AA12-AA21-AA22)/((DT)**2))-(R11+R22)/(DT)
B12 =((AA11+AA21)/(DT**2)+(R11/(DT)))
B13 =((AA12+AA22)/(DT**2)+R22/(DT))
B21 =((AA11+AA12)/(DT**2)+(R11/(DT)))
B22 =((PI1*Q1-AA11)/(DT**2)-(R11/(DT)))
B23 =(-AA12/(DT**2))
B31 =((AA21+AA22)/(DT**2)+(R22/(DT)))
B32 =(-AA21/(DT**2))
B33 =((PI2*Q2-AA22)/(DT**2)-(R22/(DT)))
write(*,*)B11,B12,B13,B21,B22,B23,B31,B32,B33
!!RHS combine B&C : inertial and stress coefficients for T+1 term!!
A(1,1)=B11*C11+B12*C21+B13*C31 !BS
A(1,2)=-Q1*(B11*C12+B12*C22+B13*C32)
A(1,3)=-Q2*(B11*C13+B12*C23+B13*C33)
A(2,1)=B21*C11+B22*C21+B23*C31 !B1
A(2,2)=-Q1*(B21*C12+B22*C22+B23*C32)
A(2,3)=-Q2*(B21*C13+B22*C23+B23*C33)
A(3,1)= B31*C11+B32*C21+B33*C31 !B2
A(3,2)=-Q1*(B31*C12+B32*C22+B33*C32)
A(3,3)=-Q2*(B31*C13+B32*C23+B33*C33)
write(*,*)A
!!setting B.C!!
DO NY=1,JM
W=1
W1=2*W*pi
W2=W1*DT*W
AA=10**8
BB=2*10**7
!left side
PS(1,NY)=AA*COS(W2)
P1(1,NY)=BB*COS(W2)
P2(1,NY)=BB*COS(W2)-4.8*10**3
END DO
!!LHS--T, T-1 time-level !!
DO NT=1,1
PS(NX,NY)=0 ! I.C:NT=1
P1(NX,NY)=0 ! I.C:NT=1
P2(NX,NY)=-4.8*10**3 !I.C:NT=1
END DO
! computing for BS !
10 DO NT=2,(OUTDT/DT)+1
BS(NT+2)=BS(NT)+BS(NT+1)
DO NX=2,KI
DO NY=2,JM
BS(NX,NY)=(C11*PS(NX,NY)-C12*Q1*P1(NX,NY)-C13*Q2*P2(NX,NY))*((2*PIS*QS-2*AA11
-2*AA12-2*AA21-2*AA22)/(DT**2)-(R11+R22)/DT)+&
(C21*PS(NX,NY)-C22*Q1*P1(NX,NY)-C23*Q2*P2(NX,NY))*((2*AA11+2*AA21)
/(DT**2)+R11/DT)+&
(C31*PS(NX,NY)-C32*Q1*P1(NX,NY)-C33*Q2*P2(NX,NY))*((2*AA12+2*AA22)
/(DT**2)+R22/DT)+&
-(C11*PS(NX,NY)-C12*Q1*P1(NX,NY)-C13*Q2*P2(NX,NY))*((PIS*QS-AA11-AA12
-AA21-AA22)/(DT**2))+&
-(C21*PS(NX,NY)-C22*Q1*P1(NX,NY)-C23*Q2*P2(NX,NY))*((AA11+AA21)
/(DT**2))+&
-(C31*PS(NX,NY)-C32*Q1*P1(NX,NY)-C33*Q2*P2(NX,NY))*((AA12+AA22)
/(DT**2))+&
(U11+4/3*G)*((C11*PS(NX+1,NY)-C12*Q1*P1(NX+1,NY)-C13*Q2*P2(NX+1,NY))-&
(2*(C11*PS(NX,NY)-C12*Q1*P1(NX,NY)-C13*Q2*P2(NX,NY)))+&
(C11*PS(NX-1,NY)-C12*Q1*P1(NX-1,NY)-C13*Q2*P2(NX-1,NY)))/(DX**2)+&
U12*((C21*PS(NX+1,NY)-C22*Q1*P1(NX+1,NY)-C23*Q2*P2(NX+1,NY))-&
(2*(C21*PS(NX,NY)-C22*Q1*P1(NX,NY)-C23*Q2*P2(NX,NY)))+&
(C21*PS(NX-1,NY)-C22*Q1*P1(NX-1,NY)-C23*Q2*P2(NX-1,NY)))/(DX**2)+&
U13*((C31*PS(NX+1,NY)-C32*Q1*P1(NX+1,NY)-C33*Q2*P2(NX+1,NY))-&
(2*(C31*PS(NX,NY)-C32*Q1*P1(NX,NY)-C33*Q2*P2(NX,NY)))+&
(C31*PS(NX-1,NY)-C32*Q1*P1(NX-1,NY)-C33*Q2*P2(NX-1,NY)))/(DX**2)+&
(U11+4/3*G)*((C11*PS(NX,NY+1)-C12*Q1*P1(NX,NY+1)-C13*Q2*P2(NX,NY+1))-&
(2*(C11*PS(NX,NY)-C12*Q1*P1(NX,NY)-C13*Q2*P2(NX,NY)))+&
(C11*PS(NX,NY-1)-C12*Q1*P1(NX,NY-1)-C13*Q2*P2(NX,NY-1)))/(DY**2)+&
U12*((C21*PS(NX,NY+1)-C22*Q1*P1(NX,NY+1)-C23*Q2*P2(NX,NY+1))-&
(2*(C21*PS(NX,NY)-C22*Q1*P1(NX,NY)-C23*Q2*P2(NX,NY)))+&
(C21*PS(NX,NY-1)-C22*Q1*P1(NX,NY-1)-C23*Q2*P2(NX,NY-1)))/(DY**2)+&
U13*((C31*PS(NX,NY+1)-C32*Q1*P1(NX,NY+1)-C33*Q2*P2(NX,NY+1))-&
(2*(C31*PS(NX,NY)-C32*Q1*P1(NX,NY)-C33*Q2*P2(NX,NY)))+&
(C31*PS(NX,NY-1)-C32*Q1*P1(NX,NY-1)-C33*Q2*P2(NX,NY-1)))/(DY**2)
BS(NT)=BS(NT-1)
BS(NT+1)=BS(NT)
END DO
END DO
! stability analysis for BS !
ERROR1=0.0
DO NT=2,(OUTDT/DT)+1
ERROR1=ERROR1+ABS(BS(NT+1)-BS(NT))
IF (ERROR1 .LT. 0.01) THEN
WRITE(1,*) "BS(NT)","T=",(T-1)*DT
ELSE
DO NT=2,(OUTDT/DT)+1
BS(NT+1)=BS(NT)
END DO
GO TO 10
END IF
END DO
END DO
! computing for BS !
12 DO NT=,(OUTDT/DT)+1
B1(NT+2)=B1(NT)+B1(NT+1)
DO NX=2,KI
DO NY=2,JM
B1(NX,NY)=(C11*PS(NX,NY)-C12*Q1*P1(NX,NY)-C13*Q2*P2(NX,NY))*((2*AA11
+2*AA12)/(DT**2)+R11/DT)+&
(C21*PS(NX,NY)-C22*Q1*P1(NX,NY)-C23*Q2*P2(NX,NY))*((2*PI1*Q1
-2*AA11)/(DT**2)-R11/DT)-&
(C31*PS(NX,NY)-C32*Q1*P1(NX,NY)-C33*Q2*P2(NX,NY))*(2*AA12
/(DT**2))-&
(C11*PS(NX,NY)-C12*Q1*P1(NX,NY)-C13*Q2*P2(NX,NY))*((AA11+AA12)
/(DT**2))-&
(C21*PS(NX,NY)-C22*Q1*P1(NX,NY)-C23*Q2*P2(NX,NY))*((PI1*Q1
-AA11)/(DT**2))-&
(C31*PS(NX,NY)-C32*Q1*P1(NX,NY)-C33*Q2*P2(NX,NY))*(-AA12
/(DT**2))+&
U21*((C11*PS(NX+1,NY)-C12*Q1*P1(NX+1,NY)-C13*Q2*P2(NX+1,NY))-&
(2*(C11*PS(NX,NY)-C12*Q1*P1(NX,NY)-C13*Q2*P2(NX,NY)))+&
(C11*PS(NX-1,NY)-C12*Q1*P1(NX-1,NY)-C13*Q2*P2(NX-1,NY)))
/(DX**2)+&
U22*((C21*PS(NX+1,NY)-C22*Q1*P1(NX+1,NY)-C23*Q2*P2(NX+1,NY))-&
(2*(C21*PS(NX,NY)-C22*Q1*P1(NX,NY)-C23*Q2*P2(NX,NY)))+&
(C21*PS(NX-1,NY)-C22*Q1*P1(NX-1,NY)-C23*Q2*P2(NX-1,NY)))
/(DX**2)+&
U23*((C31*PS(NX+1,NY)-C32*Q1*P1(NX+1,NY)-C33*Q2*P2(NX+1,NY))-&
(2*(C31*PS(NX,NY)-C32*Q1*P1(NX,NY)-C33*Q2*P2(NX,NY)))+&
(C31*PS(NX-1,NY)-C32*Q1*P1(NX-1,NY)-C33*Q2*P2(NX-1,NY)))
/(DX**2)+&
U21*((C11*PS(NX,NY+1)-C12*Q1*P1(NX,NY+1)-C13*Q2*P2(NX,NY+1))-&
(2*(C11*PS(NX,NY)-C12*Q1*P1(NX,NY)-C13*Q2*P2(NX,NY)))+&
(C11*PS(NX,NY-1)-C12*Q1*P1(NX,NY-1)-C13*Q2*P2(NX,NY-1)))
/(DY**2)+&
U22*((C21*PS(NX,NY+1)-C22*Q1*P1(NX,NY+1)-C23*Q2*P2(NX,NY+1))-&
(2*(C21*PS(NX,NY)-C22*Q1*P1(NX,NY)-C23*Q2*P2(NX,NY)))+&
(C21*PS(NX,NY-1)-C22*Q1*P1(NX,NY-1)-C23*Q2*P2(NX,NY-1)))
/(DY**2)+&
U23*((C31*PS(NX,NY+1)-C32*Q1*P1(NX,NY+1)-C33*Q2*P2(NX,NY+1))-&
(2*(C31*PS(NX,NY)-C32*Q1*P1(NX,NY)-C33*Q2*P2(NX,NY)))+&
(C31*PS(NX,NY-1)-C32*Q1*P1(NX,NY-1)-C33*Q2*P2(NX,NY-1)))/(DY**2)
B1(NT)=B1(NT-1)
B1(NT+1)=B1(NT)
END DO
END DO
! stability analysis for B1 !
ERROR2=0.0
DO NT=2,(OUTDT/DT)+1
ERROR2=ERROR2+ABS(B1(NT+1)-B1(NT))
IF (ERROR2 .LT. 0.01) THEN
WRITE(2,*)"B1(NT)","T=",(T-1)*DT
ELSE
DO NT=2,(OUTDT/DT)+1
B1(NT+1)=B1(NT)
END DO
GO TO 12
END IF
END DO
END DO
! computing for B2 !
14 DO NT=1,(OUTDT/DT)+1
B2(NT+2)=B2(NT)+B2(NT+1)
DO NX=2,KI
DO NY=2,JM
B2(NX,NY)=(C11*PS(NX,NY)-C12*Q1*P1(NX,NY)-C13*Q2*P2(NX,NY))
*((2*AA21+2*AA22)/(DT**2)+R22/DT)-&
(C21*PS(NX,NY)-C22*Q1*P1(NX,NY)-C23*Q2*P2(NX,NY)*((2*AA21)
/(DT**2))+&
(C31*PS(NX,NY)-C32*Q1*P1(NX,NY)-C33*Q2*P2(NX,NY))*((2*PI2*Q2
-2*AA22)/(DT**2)-R22/DT)-&
(C11*PS(NX,NY)-C12*Q1*P1(NX,NY)-C13*Q2*P2(NX,NY))*((AA21+AA22)
/(DT**2))+&
(C21*PS(NX,NY)-C22*Q1*P1(NX,NY)-C23*Q2*P2(NX,NY))*(AA21/(DT**2))-&
(C31*PS(NX,NY)-C32*Q1*P1(NX,NY)-C33*Q2*P2(NX,NY))*((PI2*Q2-AA22)
/(DT**2))+&
U31*((C11*PS(NX+1,NY)-C12*Q1*P1(NX+1,NY)-C13*Q2*P2(NX+1,NY))-&
(2*(C11*PS(NX,NY)-C12*Q1*P1(NX,NY)-C13*Q2*P2(NX,NY)))+&
(C11*PS(NX-1,NY)-C12*Q1*P1(NX-1,NY)-C13*Q2*P2(NX-1,NY)))/(DX**2)+&
U32*((C21*PS(NX+1,NY)-C22*Q1*P1(NX+1,NY)-C23*Q2*P2(NX+1,NY))-&
(2*(C21*PS(NX,NY)-C22*Q1*P1(NX,NY)-C23*Q2*P2(NX,NY)))+&
(C21*PS(NX-1,NY)-C22*Q1*P1(NX-1,NY)-C23*Q2*P2(NX-1,NY)))/(DX**2)+&
U33*((C31*PS(NX+1,NY)-C32*Q1*P1(NX+1,NY)-C33*Q2*P2(NX+1,NY))-&
(2*(C31*PS(NX,NY)-C32*Q1*P1(NX,NY)-C33*Q2*P2(NX,NY)))+&
(C31*PS(NX-1,NY)-C32*Q1*P1(NX-1,NY)-C33*Q2*P2(NX-1,NY)))/(DX**2)+&
U31*((C11*PS(NX,NY+1)-C12*Q1*P1(NX,NY+1)-C13*Q2*P2(NX,NY+1))-&
(2*(C11*PS(NX,NY)-C12*Q1*P1(NX,NY)-C13*Q2*P2(NX,NY)))+&
(C11*PS(NX,NY-1)-C12*Q1*P1(NX,NY-1)-C13*Q2*P2(NX,NY-1)))/(DY**2)+&
U32*((C21*PS(NX,NY+1)-C22*Q1*P1(NX,NY+1)-C23*Q2*P2(NX,NY+1))-&
(2*(C21*PS(NX,NY)-C22*Q1*P1(NX,NY)-C23*Q2*P2(NX,NY)))+&
(C21*PS(NX,NY-1)-C22*Q1*P1(NX,NY-1)-C23*Q2*P2(NX,NY-1)))/(DY**2)+&
U33*((C31*PS(NX,NY+1)-C32*Q1*P1(NX,NY+1)-C33*Q2*P2(NX,NY+1))-&
(2*(C31*PS(NX,NY)-C32*Q1*P1(NX,NY)-C33*Q2*P2(NX,NY)))+&
(C31*PS(NX,NY-1)-C32*Q1*P1(NX,NY-1)-C33*Q2*P2(NX,NY-1)))/(DY**2)
B2(NT)=B2(NT-1)
B2(NT+1)=B2(NT)
END DO
END DO
! stability analysis for B2 !
ERROR3=0.0
DO NX=2,KI
DO NY=2,JM
ERROR3=ERROR3+ABS(B2(NT+1)-B2(NT))
IF (ERROR3 .LT. 0.01) THEN
WRITE(3,*) "B2(NT)","T=",(T-1)*DT
ELSE
DO NT=2,(OUTDT/DT)+1
B2(NT+1)=B2(NT)
END DO
GO TO 14
END IF
END DO
!invert A *B!
DO NT=2,(OUTDT/DT)+1
B=(/BS(NX,NY),B1(NX,NY),B2(NX,NY)/)
X=A .IX. B
write(4,*)"X(NT)","T=",(T-1)*DT
END DO
STOP
END
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.116.77.176
※ 編輯: Gill0920 來自: 140.116.77.176 (02/28 22:08)
1F：推 laa7352:BS宣告為二維陣列，你後面把它當一維陣列來用 03/01 16:33
2F：→ Gill0920:恩恩~就是時間的迴圈我想另外寫~那要怎麼表示呢?! 03/01 17:14
3F：→ charlesdc:先用implicit none吧 你有些變數都沒宣告 用這來提醒自 03/01 22:03
4F：→ charlesdc:已變數的型態是否正確 03/01 22:03
5F：→ charlesdc:1.應該是浮點數運算誤差 解法:有時改成real*4會成功 03/01 22:06
6F：→ charlesdc:至於為什麼我也不知道 2.BS(a.b)至於要擺啥只有你自己懂 03/01 22:07
7F：→ charlesdc:3.檢驗數可以都列出來檢核最前面跟最後面幾筆 03/01 22:10
8F：→ charlesdc:4.看不懂 剩下給熱心大解 03/01 22:11
9F：→ Gill0920:相當感謝!!!! 繼續研究!! 03/02 00:01