Jivs

03-08-2015, 12:31 PM

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