[問題]麻煩幫我看看~~非常感謝!!!

看板Fortran作者Gill0920 (不斷邁進!!)時間14年前 (2010/02/28 22:07)推噓1(1推 0噓 8→)

留言9則, 3人參與討論串1/1

這是我的程式前面算係數的部分都確認過了最後的三個大廻圈不知道要怎麼表示跟輸出以下是我的問題 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)