I'm a CFD student that is trying to develop a program where I use the Finite Volume Method to calculate the temperature distribution. I already understood the philosophy and math of this method. I guess that my problem is dealing with the Source Term.First of all the equation that I'm working with is:$$\frac{\partial }{\partial x}\left ( \Gamma \frac{\partial \phi}{\partial x} \right )+\frac{\partial }{\partial y}\left ( \Gamma \frac{\partial \phi}{\partial y} \right )+S=0$$and the result after the discretization process is:$$a_P\phi_P-a_S\phi_S-a_N\phi_N-a_E\phi_E-a_W\phi_W=b$$where:$$a_E=\frac{\Gamma e}{\Delta x\delta x_e}$$$$a_W=\frac{\Gamma w}{\Delta x\delta x_w}$$$$a_N=\frac{\Gamma n}{\Delta x\delta x_n}$$$$a_S=\frac{\Gamma s}{\Delta x\delta x_s}$$$$a_P=a_E+a_W+a_N+a_S$$$$b=S_P$$Note: $\int_{x_w}^{x_e}\int_{y_s}^{y_n}Sdydx\approx S_P \Delta x \Delta y$The coefficients $a_P,a_S,a_N,a_E,a_W$ will be the same during the iteration process, so the Source Term (that must go to zero) will only change with the update of $\phi 's$. To update $\phi 's$ I will be using the TDMA method, but something is really wrong in my program.ProgramThis is the subroutine CALCT where the coefficients are calculated and the source term must be calculated. The MODT is the subroutine that keeps inserting the boundary conditions.SUBROUTINE CALCTINCLUDE 'COMMON.INC'DO I=2,NIM1DO J=2,NJM1    AN(I,J)=K(I,J+1)/(SNS(J)*DYNP(I))    AS(I,J)=K(I,J-1)/(SNS(J)*DYPS(I))    AE(I,J)=K(I+1,J)/(SEW(J)*DXEP(I))    AW(I,J)=K(I-1,J)/(SEW(J)*DXPW(I))END DOEND DOCALL MODTRESORT=0.0DO I=2,NIM1DO J=2,NJM1    AP(I,J)=AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)    RESORT=RESORT-AN(I,J)*T(I,J+1)-AS(I,J)*T(I,J-1)-AE(I,J)*T(I+1,J)-AW(I,J)*T(I-1,J)+AP(I,J)*T(I,J)END DOEND DONSWPT=5DO  N=1,NSWPT    CALL TDMA(2,2,NIM1,NJM1,IT,JT,T)END DORETURNEND SUBROUTINETDMA ALGORITHM  SUBROUTINE TDMA(ISTART,JSTART,NI,NJ,IT,JT,PHI)  IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N)  PARAMETER (NN=100)  DIMENSION PHI(IT,JT),A(NN),B(NN),C(NN),D(NN)  COMMON/COEF/AP(NN,NN),AN(NN,NN),AS(NN,NN),AE(NN,NN),AW(NN,NN),SU(NN,NN),SP(NN,NN)  WRITE(*,*)'TDMA ACTIVATED.'  NIM1=NI-1  NJM1=NJ-1  JSTM1=JSTART-1  A(JSTM1)=0.0  !C-----COMMENCE W-E SWEEP  DO 100 I=ISTART,NIM1  C(JSTM1)=PHI(I,JSTM1)  !C-----COMMENCE S-N TRAVERSE  DO 101 J=JSTART,NJM1  !C-----ASSEMBLE TDMA COEFFICIENTS  A(J)=AN(I,J)  B(J)=AS(I,J)  C(J)=AE(I,J)*PHI(I+1,J)+AW(I,J)*PHI(I-1,J)+SU(I,J)  D(J)=AP(I,J)  !C-----CALCULATE COEFFICIENTS OF RECURRENCE FORMULA  TERM=1./(D(J)-B(J)*A(J-1))  A(J)=A(J)*TERM  101 C(J)=(C(J)+B(J)*C(J-1))*TERM  !C-----OBTAIN NEW PHI@S  DO 102 JJ=JSTART,NJM1  J=NJ+JSTM1-JJ  102 PHI(I,J)=A(J)*PHI(I,J+1)+C(J)  100 CONTINUE  RETURN  END
Programming a FVM in diffusion equation
fluid dynamics;fortran;diffusion;heat transfer