Fortran Source Code & Offline Testing

FTCS Heat Equation Stability

Core Numerical Engine in Fortran 90 • 21 total downloads

heat_equation_stability.f90

! =========================================================================
! Source File: heat_equation_stability.f90
! =========================================================================

program heat_equation_stability
    implicit none
    double precision :: alpha_in,dx,dy,dz,dt_in,k_cond,rho,cp_val,T_init,T_bc
    integer :: scheme,ios,ndim,i,npts
    double precision :: alpha,inv_dx2,Fo_max,rx,ry,rz,Fo,dt_use,dt_max
    double precision :: t_cell,t_diff,pen_depth,margin
    character(len=30) :: sch_name,stab_type,status
    character(len=10) :: order_str
    double precision :: tc,fc,rxc,pc,nc,Fo_ftcs

    read(*,*,iostat=ios) alpha_in; if(ios/=0) stop 'ERROR: alpha'
    read(*,*,iostat=ios) dx;       if(ios/=0) stop 'ERROR: dx'
    read(*,*,iostat=ios) dy;       if(ios/=0) stop 'ERROR: dy'
    read(*,*,iostat=ios) dz;       if(ios/=0) stop 'ERROR: dz'
    read(*,*,iostat=ios) dt_in;    if(ios/=0) stop 'ERROR: dt'
    read(*,*,iostat=ios) k_cond;   if(ios/=0) stop 'ERROR: k'
    read(*,*,iostat=ios) rho;      if(ios/=0) stop 'ERROR: rho'
    read(*,*,iostat=ios) cp_val;   if(ios/=0) stop 'ERROR: cp'
    read(*,*,iostat=ios) T_init;   if(ios/=0) stop 'ERROR: T_init'
    read(*,*,iostat=ios) T_bc;     if(ios/=0) stop 'ERROR: T_bc'
    read(*,*,iostat=ios) scheme;   if(ios/=0) stop 'ERROR: scheme'
    if(dx<=0d0) stop 'ERROR: dx must be positive'

    if(alpha_in>0d0)then
        alpha=alpha_in
    else
        if(k_cond<=0d0.or.rho<=0d0.or.cp_val<=0d0) stop 'ERROR: need alpha or k,rho,cp'
        alpha=k_cond/(rho*cp_val)
    end if

    ndim=1
    if(dy>0d0) ndim=2
    if(dy>0d0.and.dz>0d0) ndim=3
    if(dy<=0d0) dy=1d30
    if(dz<=0d0) dz=1d30

    inv_dx2=1d0/dx**2
    if(ndim>=2) inv_dx2=inv_dx2+1d0/dy**2
    if(ndim>=3) inv_dx2=inv_dx2+1d0/dz**2

    select case(scheme)
    case(1)
        sch_name='FTCS Explicit'
        stab_type='Conditionally Stable'
        Fo_max=1d0/(2d0*dble(ndim))
        order_str='O(dt,dx2)'
    case(2)
        sch_name='Implicit BE'
        stab_type='Unconditionally Stable'
        Fo_max=1d10
        order_str='O(dt,dx2)'
    case(3)
        sch_name='Crank-Nicolson'
        stab_type='Unconditionally Stable'
        Fo_max=1d10
        order_str='O(dt2,dx2)'
    case(4)
        sch_name='DuFort-Frankel'
        stab_type='Unconditionally Stable'
        Fo_max=1d10
        order_str='O(dt2,dx2)'
    case default
        sch_name='FTCS Explicit'
        stab_type='Conditionally Stable'
        Fo_max=1d0/(2d0*dble(ndim))
        order_str='O(dt,dx2)'
    end select

    dt_max=Fo_max/(alpha*inv_dx2)
    if(dt_max>1d20) dt_max=1d10
    if(dt_in>0d0)then
        dt_use=dt_in
    else
        dt_use=min(dt_max*0.9d0,1d5)
    end if

    rx=alpha*dt_use/dx**2; ry=0d0; rz=0d0
    if(ndim>=2) ry=alpha*dt_use/dy**2
    if(ndim>=3) rz=alpha*dt_use/dz**2
    Fo=rx+ry+rz

    Fo_ftcs=1d0/(2d0*dble(ndim))
    if(scheme>=2)then
        status='UNCONDITIONAL'; margin=100d0
    else if(Fo<=Fo_max)then
        status='STABLE'; margin=(Fo_max-Fo)/max(Fo_max,1d-30)*100d0
    else
        status='UNSTABLE'; margin=(Fo_max-Fo)/max(Fo_max,1d-30)*100d0
    end if

    t_cell=dx**2/alpha; pen_depth=2d0*sqrt(alpha*dt_use); t_diff=(10d0*dx)**2/alpha

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   HEAT EQUATION STABILITY CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUT CONDITIONS ----------------------------------------'
    write(*,'(A,ES14.6,A)') '  alpha (diffusivity)     = ', alpha, ' m2/s'
    write(*,'(A,ES14.6,A)') '  k (conductivity)        = ', k_cond, ' W/m.K'
    write(*,'(A,F12.4,A)')  '  rho                     = ', rho, ' kg/m3'
    write(*,'(A,F12.4,A)')  '  cp                      = ', cp_val, ' J/kg.K'
    write(*,'(A,F12.2,A)')  '  T_init                  = ', T_init, ' K'
    write(*,'(A,F12.2,A)')  '  T_bc                    = ', T_bc, ' K'
    write(*,*)
    write(*,'(A)') '--- GRID INFORMATION ----------------------------------------'
    write(*,'(A,ES14.6,A)') '  dx                      = ', dx, ' m'
    if(ndim>=2) write(*,'(A,ES14.6,A)') '  dy                      = ', dy, ' m'
    if(ndim>=3) write(*,'(A,ES14.6,A)') '  dz                      = ', dz, ' m'
    write(*,'(A,I2,A)') '  Dimensions              = ', ndim, 'D'
    write(*,*)
    write(*,'(A)') '--- SCHEME INFORMATION --------------------------------------'
    write(*,'(A,A)') '  Scheme                  = ', trim(sch_name)
    write(*,'(A,A)') '  Stability type          = ', trim(stab_type)
    if(Fo_max<1d5) write(*,'(A,F12.6)') '  Fo_max (limit)          = ', Fo_max
    if(Fo_max>=1d5) write(*,'(A)') '  Fo_max (limit)          = unlimited'
    write(*,'(A,A)') '  Accuracy order          = ', trim(order_str)
    write(*,*)
    write(*,'(A)') '--- TIME STEP -----------------------------------------------'
    write(*,'(A,ES14.6,A)') '  dt (used)               = ', dt_use, ' s'
    write(*,'(A,ES14.6,A)') '  dt_max (FTCS stable)    = ', min(dt_max,1d10), ' s'
    write(*,*)
    write(*,'(A)') '--- FOURIER NUMBERS -----------------------------------------'
    write(*,'(A,F14.6)') '  Fo (total)              = ', Fo
    write(*,'(A,F14.6)') '  rx = alpha*dt/dx2       = ', rx
    if(ndim>=2) write(*,'(A,F14.6)') '  ry = alpha*dt/dy2       = ', ry
    if(ndim>=3) write(*,'(A,F14.6)') '  rz = alpha*dt/dz2       = ', rz
    if(Fo_max<1d5) write(*,'(A,F14.6)') '  Fo_max (scheme)         = ', Fo_max
    write(*,*)
    write(*,'(A)') '--- STABILITY STATUS ----------------------------------------'
    write(*,'(A,A)') '  STATUS                  = ', trim(status)
    write(*,'(A,F12.2,A)') '  Margin                  = ', margin, ' %'
    write(*,*)
    write(*,'(A)') '--- DIFFUSION CHARACTERISTICS -------------------------------'
    write(*,'(A,ES14.6,A)') '  t_cell = dx2/alpha      = ', t_cell, ' s'
    write(*,'(A,ES14.6,A)') '  Penetration depth       = ', pen_depth, ' m'
    write(*,'(A,F14.4)') '  Pen. depth / dx         = ', pen_depth/dx
    write(*,'(A,ES14.6,A)') '  t_diff (10*dx domain)   = ', t_diff, ' s'
    if(dt_use>0d0) write(*,'(A,F14.1)') '  Steps to t_diff         = ', t_diff/dt_use
    write(*,*)
    write(*,'(A)') '--- ALL SCHEMES COMPARISON ----------------------------------'
    write(*,'(A)') '  Scheme              Fo_max      Stable?  Order'
    write(*,'(A)') '  --------------------------------------------------------'
    if(Fo<=Fo_ftcs)then
      write(*,'(A,F10.4,A)') '  FTCS               ',Fo_ftcs,'      YES      O(dt,dx2)'
    else
      write(*,'(A,F10.4,A)') '  FTCS               ',Fo_ftcs,'      NO       O(dt,dx2)'
    end if
    write(*,'(A)') '  Implicit BE         unlimited    YES      O(dt,dx2)'
    write(*,'(A)') '  Crank-Nicolson      unlimited    YES      O(dt2,dx2)'
    write(*,'(A)') '  DuFort-Frankel      unlimited    YES*     O(dt2,dx2)'
    write(*,*)

    write(*,'(A)') '--- PROFILE vs dt ------------------------------------------'
    write(*,'(A)') '  dt          Fo          rx          Status      pen_depth   N_steps'
    write(*,'(A)') '  --------------------------------------------------------------------------'
    npts=40
    do i=1,npts
      if(dt_max<1d5)then
        tc=dt_max/20d0+(dt_max*4d0-dt_max/20d0)*dble(i-1)/dble(npts-1)
      else
        tc=dt_use/20d0+(dt_use*40d0-dt_use/20d0)*dble(i-1)/dble(npts-1)
      end if
      if(tc<=0d0) tc=1d-15
      fc=alpha*tc*inv_dx2; rxc=alpha*tc/dx**2; pc=2d0*sqrt(alpha*tc)
      nc=t_diff/max(tc,1d-30)
      if(fc<=Fo_ftcs)then
        write(*,'(ES12.4,2X,F10.6,2X,F10.6,2X,A8,2X,ES12.4,2X,ES12.4)') tc,fc,rxc,'STABLE  ',pc,nc
      else
        write(*,'(ES12.4,2X,F10.6,2X,F10.6,2X,A8,2X,ES12.4,2X,ES12.4)') tc,fc,rxc,'UNSTABLE',pc,nc
      end if
    end do

    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') '  Fo = alpha*dt/dx^2  (Fourier number, Eq. 2.80)'
    write(*,'(A)') '  FTCS: Fo <= 1/(2*ndim)'
    write(*,'(A)') '  1D: Fo<=0.5, 2D: Fo<=0.25, 3D: Fo<=1/6'
    write(*,'(A)') '  dt_max = Fo_max / (alpha * sum(1/dxi^2))'
    write(*,'(A)') '  Penetration depth = 2*sqrt(alpha*t)'
    write(*,'(A)') '============================================================'

end program heat_equation_stability

Solver Description

Analyze numerical stability limits for 1D/2D transient diffusion solvers using the FTCS scheme.

Key Numerical Methods & Architecture

Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.

🛠️ Local Compilation

To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).

Compilation Command:

gfortran -O3 heat_equation_stability.f90 -o heat_equation_stability

Execution Command:

Execute the program by feeding the sample input file into the program using stdin redirection:

heat_equation_stability < input.txt

📥 Downloads & Local Files

Preview of the required input file (input.txt):

! a_i\nx [m]\ny [m]\nz [m]\nt [s]\nk_i\nr_i\ncp [J/kgK]\nti_i\ntb_i\nsc_i
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
! Parameter 6
0.0
! Parameter 7
0.0
! Parameter 8
0.0
! Parameter 9
0.0
! Parameter 10
0.0
! Parameter 11
8

📦 Download Offline Bundle (ZIP) 7 dl 📄 Download input.txt 8 dl 📥 Download heat_equation_stability.f90 6 dl

💻 Fortran Source Code Library

FTCS Heat Equation Stability

Solver Description

Key Numerical Methods & Architecture

🛠️ Local Compilation

📥 Downloads & Local Files