Hi All,
I am trying to figure out how to allocate the size for Array dynamically, as in I want it to be from 0 to a Number till grid independent solution is reached by a Condition.
So below is my program, I run it many time for different values of N. But then I don't want to do this, could you please help me with this
Thank You!
P.S: I have tried different ways but then since the entire program depends on N. I am not able to figure out what to do!
program HW1_P1_a
Integer, Parameter :: N=250 !Globally used variable, hence PARAMETER
Double Precision Del_X,A(N),B(N),C(N),D(N),PA,QA,RA,PB,QB,RB
Double Precision L,M,Theta_B,Theta_L, THETA(N)
Integer I,T_B,T_L,T_INF
! "Details about Arguments used: All must be Double Precision (except N)"
! N Number of mesh points
! Del_X Mesh size
! A() Lower diagonal of tridiagonal matrix
! B() Main diagonal of tridiagonal matrix
! C() Upper diagonal of tridiagonal matrix
! D() Right side vector (do not adjust for bcs)
! PA Left b.c.(K=1): Coefficient of Theta
! QA Left b.c.(K=1): Coefficient of dTheta/Del_X
! RA Left b.c.(K=1): Right hand side
! PB Right b.c.(K=N): Coefficient of Theta
! QB Right b.c.(K=N): Coefficient of dTheta/Del_X
! RB Right b.c.(K=N): Right hand side
! L Length of Fin
! T_B Temperature at Base of Fin
! T_L Temperature at the Tip of Fin
! T_inf Temperature of surrounding air
! Theta_B Temperature excess at Fin Base
! Theta_L Temperature excess at Fin Tip
! P Perimeter
! h Heat Transfer Coeffecient
! Ac Cross Sectional Area of Fin
! K Thermal Conductivity
! Inititalising values of the variables, "Given Conditions"
T_B=200
T_L=60
T_INF=20
L=0.1
P=0.09
Ac=7.E-4
h=15
K=80
! Equation for Excess Temperature
Theta_B=T_B-T_INF
Theta_L=T_L-T_INF
! Equation for Delta X and M
Del_X=L/(N-1) !Mesh size
M=SQRT((h*P)/(K*Ac))
! Fin equation derivatives are approximated using second order accurate finite difference and multiplying with (Del_X)^2
g=M**2
f=0
! Elements of the Matrix
DO I=1,N
A(i) = 1 - (Del_X/2)*f
B(i) = -(2 + (Del_X**2)*g)
C(i) = 1 + (Del_X/2)*f
D(i) = 0
ENDDO
! Dirchlet Boundary Condition at the Base
PA=1
QA=0
RA=Theta_B
! Dirclet Boundary condition at the Tip
PB=1
QB=0
RB=Theta_L
call THOMAS(N,Del_X,A,B,C,D,PA,QA,RA,PB,QB,RB,THETA)
Write(*,*) 'Values of Theta by Numerical Solutions for N:',N
DO I=1,N
write (*,*) (I*Del_X),THETA(I)
ENDDO
End HW1_P1_a
!**********************************************************************
!
! SUBROUTINE: THOMAS
!
! DESCRIPTION: Use the Thomas algorithm to solve a system of
! linear algebraic equations where the
! coefficient matrix is tridiagonal. Boundary
! conditions can be Dirichlet, Neumann or Mixed.
!
! Left boundary condition is (X = Xmin):
!
! PA*Theta + QA*dTheta/Del_X = RA
!
! Right boundary condition is (X = Xmax):
!
! PB*Theta + QB*dTheta/Del_X = RB
!
! System of Equations:
!
! B1 C1 0 0 . . 0 = D1
! A2 B2 C2 0 . .
! 0 A3 B3 C3 . .
! . . . .
! . Ai Bi Ci = Di
! 0 . . . . An Bn = Dn
!
!
! RETURN ARGUMENTS:
! THETA() Function value
!
! VARIABLES:
!
! **********************************************************************
!
SUBROUTINE THOMAS(N,Del_X,A,B,C,D,PA,QA,RA,PB,QB,RB,THETA)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (KDIM=10001)
DIMENSION A(*),B(*),C(*),D(*),THETA(*)
DIMENSION F(KDIM),DELTA(KDIM)
! -- Adjust coefficients for boundary conditions
B(1) = QA*B(1) + 2.D0*Del_X*PA*A(1)
C(1) = QA*(C(1) + A(1))
D(1) = QA*D(1) + 2.D0*Del_X*RA*A(1)
A(N) = QB*(A(N) + C(N))
B(N) = QB*B(N) - 2.D0*Del_X*PB*C(N)
D(N) = QB*D(N) - 2.D0*Del_X*RB*C(N)
! -- Forward elimination
F(1) = C(1)/B(1)
DELTA(1) = D(1)/B(1)
DO K=2,N
X1 = B(K) - A(K)*F(K-1)
X1 = 1.D0/X1
F(K) = C(K)*X1
DELTA(K) = X1*(D(K) - A(K)*DELTA(K-1))
END DO
! -- Back substitution
THETA(N) = DELTA(N)
DO K=N-1,1,-1
THETA(K) = DELTA(K) - F(K)*THETA(K+1)
END DO
RETURN
END
I am trying to figure out how to allocate the size for Array dynamically, as in I want it to be from 0 to a Number till grid independent solution is reached by a Condition.
So below is my program, I run it many time for different values of N. But then I don't want to do this, could you please help me with this
Thank You!
P.S: I have tried different ways but then since the entire program depends on N. I am not able to figure out what to do!
program HW1_P1_a
Integer, Parameter :: N=250 !Globally used variable, hence PARAMETER
Double Precision Del_X,A(N),B(N),C(N),D(N),PA,QA,RA,PB,QB,RB
Double Precision L,M,Theta_B,Theta_L, THETA(N)
Integer I,T_B,T_L,T_INF
! "Details about Arguments used: All must be Double Precision (except N)"
! N Number of mesh points
! Del_X Mesh size
! A() Lower diagonal of tridiagonal matrix
! B() Main diagonal of tridiagonal matrix
! C() Upper diagonal of tridiagonal matrix
! D() Right side vector (do not adjust for bcs)
! PA Left b.c.(K=1): Coefficient of Theta
! QA Left b.c.(K=1): Coefficient of dTheta/Del_X
! RA Left b.c.(K=1): Right hand side
! PB Right b.c.(K=N): Coefficient of Theta
! QB Right b.c.(K=N): Coefficient of dTheta/Del_X
! RB Right b.c.(K=N): Right hand side
! L Length of Fin
! T_B Temperature at Base of Fin
! T_L Temperature at the Tip of Fin
! T_inf Temperature of surrounding air
! Theta_B Temperature excess at Fin Base
! Theta_L Temperature excess at Fin Tip
! P Perimeter
! h Heat Transfer Coeffecient
! Ac Cross Sectional Area of Fin
! K Thermal Conductivity
! Inititalising values of the variables, "Given Conditions"
T_B=200
T_L=60
T_INF=20
L=0.1
P=0.09
Ac=7.E-4
h=15
K=80
! Equation for Excess Temperature
Theta_B=T_B-T_INF
Theta_L=T_L-T_INF
! Equation for Delta X and M
Del_X=L/(N-1) !Mesh size
M=SQRT((h*P)/(K*Ac))
! Fin equation derivatives are approximated using second order accurate finite difference and multiplying with (Del_X)^2
g=M**2
f=0
! Elements of the Matrix
DO I=1,N
A(i) = 1 - (Del_X/2)*f
B(i) = -(2 + (Del_X**2)*g)
C(i) = 1 + (Del_X/2)*f
D(i) = 0
ENDDO
! Dirchlet Boundary Condition at the Base
PA=1
QA=0
RA=Theta_B
! Dirclet Boundary condition at the Tip
PB=1
QB=0
RB=Theta_L
call THOMAS(N,Del_X,A,B,C,D,PA,QA,RA,PB,QB,RB,THETA)
Write(*,*) 'Values of Theta by Numerical Solutions for N:',N
DO I=1,N
write (*,*) (I*Del_X),THETA(I)
ENDDO
End HW1_P1_a
!**********************************************************************
!
! SUBROUTINE: THOMAS
!
! DESCRIPTION: Use the Thomas algorithm to solve a system of
! linear algebraic equations where the
! coefficient matrix is tridiagonal. Boundary
! conditions can be Dirichlet, Neumann or Mixed.
!
! Left boundary condition is (X = Xmin):
!
! PA*Theta + QA*dTheta/Del_X = RA
!
! Right boundary condition is (X = Xmax):
!
! PB*Theta + QB*dTheta/Del_X = RB
!
! System of Equations:
!
! B1 C1 0 0 . . 0 = D1
! A2 B2 C2 0 . .
! 0 A3 B3 C3 . .
! . . . .
! . Ai Bi Ci = Di
! 0 . . . . An Bn = Dn
!
!
! RETURN ARGUMENTS:
! THETA() Function value
!
! VARIABLES:
!
! **********************************************************************
!
SUBROUTINE THOMAS(N,Del_X,A,B,C,D,PA,QA,RA,PB,QB,RB,THETA)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (KDIM=10001)
DIMENSION A(*),B(*),C(*),D(*),THETA(*)
DIMENSION F(KDIM),DELTA(KDIM)
! -- Adjust coefficients for boundary conditions
B(1) = QA*B(1) + 2.D0*Del_X*PA*A(1)
C(1) = QA*(C(1) + A(1))
D(1) = QA*D(1) + 2.D0*Del_X*RA*A(1)
A(N) = QB*(A(N) + C(N))
B(N) = QB*B(N) - 2.D0*Del_X*PB*C(N)
D(N) = QB*D(N) - 2.D0*Del_X*RB*C(N)
! -- Forward elimination
F(1) = C(1)/B(1)
DELTA(1) = D(1)/B(1)
DO K=2,N
X1 = B(K) - A(K)*F(K-1)
X1 = 1.D0/X1
F(K) = C(K)*X1
DELTA(K) = X1*(D(K) - A(K)*DELTA(K-1))
END DO
! -- Back substitution
THETA(N) = DELTA(N)
DO K=N-1,1,-1
THETA(K) = DELTA(K) - F(K)*THETA(K+1)
END DO
RETURN
END
