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