fdm0308.f90

2d FDM 源代码 用于均匀磁场的位场延拓

PROGRAM FDM 
!***************************************************************!
!说明：U0为磁场初值  A为系数矩阵 B为右端项 JS 解方程组中的工作数组  
!共六行 九列  总节点数：6*9=54
IMPLICIT NONE
INTEGER I,L
REAL U0(6,9),A(54,54),B(54),C(6,9,54)
REAL JS(54)             
CALL READIN(I,U0)         !读入第一边界初始值
CALL INT_B(U0,B)          !初始化B值
CALL INITI_C(C)           !初始化C矩阵 
CALL CACUL_C(C)           !计算C系数矩阵
CALL WRITE_A(C)
CALL TRANSFORM_A(C,A)     !将C系数矩阵转换为二维矩阵A(54,54)     
CALL AGGJE(A,54,B,L,JS)   !用全选主元的高斯-亚当消去法接线性方程组 详见《徐世良算法集》P5!
IF (L.NE.0) THEN		 
	  OPEN(3,FILE='OUT.TXT',POSITION='APPEND')
	  DO I=1,54
	    WRITE(3,10)I,B(I)
	    WRITE(*,10)I,B(I)
	  END DO
	END IF
	CLOSE(3)
10	FORMAT(1X,'U',I4,'=',F8.3)     !将计算结果输出
END




!从文件中读入第一边界(地面)初始值
SUBROUTINE READIN(I,U0)
REAL U0(6,9)
OPEN(1,FILE='INDATAU0.TXT')
READ(1,*)
READ(1,*)
READ(1,*)(U0(1,I),I=1,9) 
READ(1,*)(U0(6,I),I=1,9)
READ(1,*)(U0(I,1),I=2,5)
READ(1,*)(U0(I,9),I=2,5)
CLOSE(1)
RETURN
END


!初始化B
SUBROUTINE INT_B(U0,B)
INTEGER I,J,K
REAL B(54),U0(6,9)
DO I=1,6
    DO J=1,9
	 K=1
     B(K)=U0(I,J)
	 K=K+1
    END DO
END DO 
RETURN
END


!初始化C矩阵
SUBROUTINE INITI_C(C)
real C(6,9,54)
INTEGER I,J,K
DO K=1,54
  DO I=1,6
    DO J=1,9
	  C(I,J,K)=0.
	END DO
  END DO
END DO
!地面边界节点初始化
C(1,1,1)=1.
C(1,2,2)=1.
C(1,3,3)=1.
C(1,4,4)=1.
C(1,5,5)=1.
C(1,6,6)=1.
C(1,7,7)=1.
C(1,8,8)=1.
C(1,9,9)=1.
!左侧边界节点初始化
I=2        
DO K=10,36,9
  C(I,1,K)=1.
  I=I+1
END DO
!右侧边界节点初始化
I=2
DO K=18,54,9
C(I,9,K)=1.
I=I+1
 END DO
!上部边界节点初始化
K=46
J=0
DO WHILE(K<=54)
  J=J+1
  K=K+1
  C(6,J,K)=1.
END DO 
RETURN
END			!所有节点C矩阵初始化完毕	   

!计算C矩阵系数 
SUBROUTINE CACUL_C(C) 
REAL C(6,9,54)
REAL H
INTEGER I,PL,PU,J
I=1
DO PL=11,38,9        !从第二行开始逐行计算C(I,I,K)
I=I+1
PU=PL+6
J=2
H=1.2
CALL CACUL_C0(C,H,I,PL,PU,J)
 END DO
RETURN
END

!计算C子矩阵
SUBROUTINE 	CACUL_C0(C,H,I,PL,PU,J)
REAL C(6,9,54)
INTEGER M,K,J,PL,PU
REAL A,B
DO K=PL,PU
      M=ABS(5-J)
	  A=1./((H**(2*M+1))*(1.+H))
	  B=1./((H**(2*I-1))*(1.+H))
  IF(J<5)THEN
      C(I,J-1,K)=A
	  C(I,J+1,K)=A*H
	  C(I-1,J,K)=B*H
	  C(I+1,J,K)=B
	  C(I,J,K)=-1.*(C(I,J-1,K)+C(I,J+1,K)+C(I-1,J,K)+C(I+1,J,K))
  ELSE IF(J>5)THEN
      C(I,J-1,K)=A*H
	  C(I,J+1,K)=A
	  C(I-1,J,K)=B*H
	  C(I+1,J,K)=B
	  C(I,J,K)=-1.*(C(I,J-1,K)+C(I,J+1,K)+C(I-1,J,K)+C(I+1,J,K))
  ELSE
	  C(I,J-1,K)=1./(2*(H**2))
	  C(I,J+1,K)=1./(2*(H**2))
	  C(I-1,J,K)=B*H
	  C(I+1,J,K)=B
	  C(I,J,K)=-1.*(C(I,J-1,K)+C(I,J+1,K)+C(I-1,J,K)+C(I+1,J,K))
  END IF
      J=J+1
END DO
RETURN
END

!把C矩阵打印在A.TXT文件上
SUBROUTINE WRITE_A(C)
REAL C(6,9,54)
OPEN(2,FILE='A.TXT')
Do K=1,54
 WRITE(2,'(2X,225F8.5)')((C(I,J,K),J=1,9),I=1,6)
END DO
CLOSE(2)
END

!把3D矩阵C(6,9,54)转换为2D矩阵A(54,54)
SUBROUTINE TRANSFORM_A(C,A)
REAL C(6,9,54)
REAL A(54,54)
INTEGER I,J,K
INTEGER N
DO K=1,54
 N=0
  DO I=1,6
    DO J=1,9
	  N=N+1
	  A(K,N)=C(I,J,K)
	END DO
   END DO
END DO
RETURN
END

!用高斯-约当消去法解线性方程组
	SUBROUTINE AGGJE(A,N,B,L,JS)
	REAL A(N,N),B(N),JS(N)
	REAL T,D
	INTEGER L,K,J,I,IS
	L=1
	DO 100 K=1,N
	  D=0.0
	  DO 10 I=K,N
	  DO 10 J=K,N
	    IF (ABS(A(I,J)).GT.D) THEN
	      D=ABS(A(I,J))
	      JS(K)=J
	      IS=I
	    END IF
10	  CONTINUE
	  IF (D+1.0.EQ.1.0) THEN
	    WRITE(*,20)
	    L=0
	    RETURN
	  END IF
20	  FORMAT(1X,'  FAIL  ')
	  DO 30 J=K,N
	    T=A(K,J)
	    A(K,J)=A(IS,J)
	    A(IS,J)=T
30	  CONTINUE
	  T=B(K)
	  B(K)=B(IS)
	  B(IS)=T
	  DO 50 I=1,N
	    T=A(I,K)
	    A(I,K)=A(I,JS(K))
	    A(I,JS(K))=T
50	  CONTINUE
	  T=A(K,K)
	  DO 60 J=K+1,N
	    IF (A(K,J).NE.0.0) A(K,J)=A(K,J)/T
60	  CONTINUE
	  B(K)=B(K)/T
	  DO 80 J=K+1,N
	    IF (A(K,J).NE.0.0) THEN
	      DO 70 I=1,N
	        IF ((I.NE.K).AND.(A(I,K).NE.0.0)) THEN
                  A(I,J)=A(I,J)-A(I,K)*A(K,J)
	        END IF
70	      CONTINUE
	    END IF
80	  CONTINUE
	  DO 90 I=1,N
	    IF ((I.NE.K).AND.(A(I,K).NE.0.0)) THEN
              B(I)=B(I)-A(I,K)*B(K)
	    END IF
90	  CONTINUE
100	CONTINUE
	DO 110 K=N,1,-1
	  IF (K.NE.JS(K)) THEN
	    T=B(K)
	    B(K)=B(JS(K))
	    B(JS(K))=T
	  END IF
110	CONTINUE
	RETURN
	END