! File: solve_heat_diffusion_v1.f90
! Purpose: OOP implementation in Fortran of the ADE method of Barakat & Clark
! (1966), for solving the time-dependent heat-diffusion equation in 2D,
! with a linearly-varying temperature profile on the domain-boundaries.
!NOTES: - (for gfortran users): gfortran-4.8 still had no support for
! 'final'-methods, and an error will be raised by that version (or
! earlier) of that compiler; HOWEVER, the final-procedure can be safely
! removed in our case, since it is included for demonstration-purposes
! only (the arrays which it de-allocates are freed anyway by the runtime
! system, when the 'Solver'-instance goes out of scope). Comment
! corresponding line below, to disable use of this feature.
! - The file "config_file_formatted.in", containing the parameters for the
! solver, needs to be located in working directory.
! - Lines containing only comments with dots can be ignored (they just
! help with including code in TeX-files).
17 module NumericKinds
18 implicit none
! KINDs for different types of REALs
21 integer, parameter :: &
R_SP = selected_real_kind( 6, 37 ), &
R_DP = selected_real_kind( 15, 307 ), &
R_QP = selected_real_kind( 33, 4931 )
! Alias for precision that we use in the program (change this to any of the
! values 'R_SP', 'R_DP', or 'R_QP', to switch to another precision globally).
27 integer, parameter :: RK = R_DP ! if changing this, also change RK_FMT
! KINDs for different types of INTEGERs
30 integer, parameter :: &
I1B = selected_int_kind(2), & ! max = 127
I2B = selected_int_kind(4), & ! max ~ 3.28x10^4
I3B = selected_int_kind(9), & ! max ~ 2.15x10^9
I4B = selected_int_kind(18) ! max ~ 9.22x10^18
! Alias for integer-precision (analogue role to RK above).
36 integer, parameter :: IK = I3B
! Edit-descriptors for real-values
39 character(len=*), parameter :: R_SP_FMT = "f0.6", &
R_DP_FMT = "f0.15", R_QP_FMT = "f0.33"
! Alias for output-precision to use in the program (keep this in sync with RK)
42 character(len=*), parameter :: RK_FMT = R_DP_FMT
43 end module NumericKinds
45 module Config_class
46 use NumericKinds
47 implicit none
48 private
50 type, public :: Config
51 real(RK) :: mDiffusivity = 1.15E-6_RK, & ! sandstone
! NOTE: "physical" units here (Celsius)
mTempA = 100._RK, &
mTempB = 75._RK, &
mTempC = 50._RK, &
mTempD = 25._RK, &
mSideLength = 30._RK
58 integer(IK) :: mNx = 200 ! # of points for square side-length
59 end type Config
! Generic IFACE for user-defined CTOR
62 interface Config
63 module procedure createConfig
64 end interface Config
66 contains
67 type(Config) function createConfig( cfgFilePath )
68 character(len=*), intent(in) :: cfgFilePath
69 integer :: cfgFileID
70 open( newunit=cfgFileID, file=trim(cfgFilePath), status='old', action='read' )
! .....................................................
72 read(cfgFileID,*) createConfig%mTempA
73 read(cfgFileID,*) createConfig%mTempB
74 read(cfgFileID,*) createConfig%mTempC
75 read(cfgFileID,*) createConfig%mTempD
76 read(cfgFileID,*) createConfig%mNx
77 read(cfgFileID,*) createConfig%mDiffusivity
78 read(cfgFileID,*) createConfig%mSideLength
79 close(cfgFileID)
80 end function createConfig
81 end module Config_class
83 module Solver_class
84 use NumericKinds
85 use Config_class
86 implicit none
87 private
89 type, public :: Solver
90 private ! Hide internal-data from users.
91 type(Config) :: mConfig
92 real(RK) :: mNt, & ! # of iterations to simulate a characteristic time
mDx, mDt, mA, mB ! Configuration-dependent factors.
94 real(RK), allocatable, dimension(:,:) :: mU, mV ! main work-arrays
95 integer(IK) :: mNumItersMax, mCurrIter = 0
96 contains
97 private ! By default, hide methods (and expose as needed).
98 procedure, public :: init
99 procedure, public :: run
100 procedure, public :: writeAscii
101 procedure, public :: getTemp
! Internal methods (users don't need to know about these).
103 procedure :: advanceU
104 procedure :: advanceV
!final :: cleanup ! NOTE: may need to comment-out for gfortran!
106 end type Solver
108 contains
subroutine init( this, cfgFilePath, simTime ) ! initialization subroutine
110 class(Solver), intent(inout) :: this
111 character(len=*), intent(in) :: cfgFilePath
112 real(RK), intent(in) :: simTime
! .....................................................
114 integer(IK) :: nX, i, j
115 real(RK) :: lambda
117 this%mConfig = Config( cfgFilePath ) ! call component CTOR
118 nX = this%mConfig%mNx ! for making code below more compact
! conservative choice for N_t, to resolve transients -- see Barakat & Clark (1966)
120 this%mNt = nX**2
! evaluate derived parameter 'mLambda' in Solver, based on configuration
122 this%mNumItersMax = int( simTime*this%mNt )
123 this%mDx = 1. / nX
124 this%mDt = 1. / this%mNt
125 lambda = (2.*this%mDt) / (this%mDx**2)
126 this%mA = (1.-lambda)/(1.+lambda)
127 this%mB = lambda/(2.*(1.+lambda))
! allocate memory for internal arrays
129 allocate( this%mU(0:nX, 0:nX), this%mV(0:nX, 0:nX) )
! initialize mU-field:
! - set initial temperature everywhere...
132 this%mU = 1.
! - BUT re-write @ boundaries, for correct BCs
! -- North
135 this%mU(:, nX) = [ (1./3.*(i/real(nX, RK))+2./3., i=0,nX) ]
! -- West
137 this%mU(0, :) = [ (1./3.*(j/real(nX, RK))+1./3., j=0,nX) ]
! -- South
139 this%mU(:, 0) = [ (-1./3.*(i/real(nX, RK))+1./3., i=0,nX) ]
! -- East
141 this%mU(nX, :) = [ (j/real(nX, RK), j=0,nX) ]
! initialize mV-field (from mU-field)
143 this%mV = this%mU
144 end subroutine init
146 real(RK) function getTemp( this, i, j ) ! GETter for temperature
147 class(Solver), intent(in) :: this
148 integer(IK), intent(in) :: i, j
149 getTemp = 0.5*( this%mU(i,j) + this%mV(i,j) )
150 end function getTemp
subroutine run( this ) ! method for time-marching
153 class(Solver), intent(inout) :: this
154 integer(IK) :: k ! dummy index (time-marching)
156 do k=1, this%mNumItersMax ! MAIN loop
! simple progress-monitor
158 if( mod(k-1, (this%mNumItersMax-1)/10) == 0 ) then
159 write(*, '(i5,a)') nint((k*100.0)/this%mNumItersMax), "%"
160 end if
! defer work to private methods
162 call this%advanceU()
163 call this%advanceV()
164 this%mCurrIter = this%mCurrIter + 1 ! tracking time step
165 end do
166 end subroutine run
168 subroutine advanceU( this )
169 class(Solver), intent(inout) :: this
170 integer(IK) :: i, j ! local variables
! actual update for 'mU'-field (NE-ward)
172 do j=1, this%mConfig%mNx-1 ! do NOT update
173 do i=1, this%mConfig%mNx-1 ! boundaries
174 this%mU(i,j) = this%mA*this%mU(i,j) + this%mB*( &
this%mU(i-1,j) + this%mU(i+1,j) + this%mU(i,j-1) + this%mU(i,j+1) )
176 end do
177 end do
178 end subroutine advanceU
180 subroutine advanceV( this ) ! similar to 'advanceU'
181 class(Solver), intent(inout) :: this
! .....................................................
183 integer(IK) :: i, j ! local variables
! actual update for 'mV'-field (SW-ward)
185 do j=this%mConfig%mNx-1, 1, -1 ! do NOT update
186 do i=this%mConfig%mNx-1, 1, -1 ! boundaries
187 this%mV(i,j) = this%mA*this%mV(i,j) + this%mB*( &
this%mV(i-1,j) + this%mV(i+1,j) + this%mV(i,j-1) + this%mV(i,j+1) )
189 end do
190 end do
191 end subroutine advanceV
! method for producing a ASCII output file
194 subroutine writeAscii( this, outFilePath )
195 class(Solver), intent(in) :: this
196 character(len=*), intent(in) :: outFilePath
! .....................................................
198 integer(IK) :: x, y, outFileID ! local variables
200 open( newunit=outFileID, file=trim(outFilePath), status='replace', action='write' )
201 write(outFileID, '(a)') &
"# output file for program solve_heat_diffusion_v1.f90"
203 write(outFileID, '(a,2x,a)') '"s"', "# time unit"
204 write(outFileID, '(f0.8,2x,a)') &
(this%mCurrIter*this%mConfig%mSideLength**2)/ &
(this%mConfig%mDiffusivity*this%mNt), &
"# current time"
208 write(outFileID, '(a,2x,a)') '"m"', "# X unit"
209 write(outFileID, '(i0,2x,a)') this%mConfig%mNx, "# Nx"
210 write(outFileID, '(a,2x,a)') '"m"', "# Y unit"
211 write(outFileID, '(i0,2x,a)') this%mConfig%mNx, "# Ny"
212 write(outFileID, '(a,2x,a)') '"degree~C"', "# temperature unit"
! X-axis
214 do x=0, this%mConfig%mNx
215 write(outFileID, '(f0.8,2x)', advance='no') this%mDx*this%mConfig%mSideLength*x
216 end do
217 write(outFileID, '(a)') "# XVals"
! Y-axis
219 do y=0, this%mConfig%mNx
220 write(outFileID, '(f0.8,2x)', advance='no') this%mDx*this%mConfig%mSideLength*y
221 end do
222 write(outFileID, '(a)') "# YVals"
! simulation results
224 write(outFileID, '(a)') "# from next line to end: simulated temperature"
225 do y=0, this%mConfig%mNx
226 do x=0, this%mConfig%mNx
227 write(outFileID, '(f0.8,2x)', advance='no') &
this%mConfig%mTempD+this%getTemp(x,y)*(this%mConfig%mTempA-this%mConfig%mTempD)
229 end do
230 write(outFileID, *) ! newline to separate rows for R visualization script
231 end do
232 close(outFileID)
233 end subroutine writeAscii
! destructor method
236 subroutine cleanup( this )
! 'class' -> 'type' (dummy-arg cannot be polymorphic for final procedures)
238 type(Solver), intent(inout) :: this
! in this version, we only deallocate memory
240 deallocate( this%mU, this%mV )
241 end subroutine cleanup
242 end module Solver_class
program solve_heat_diffusion_v1
use NumericKinds
use Solver_class
247 implicit none
type(Solver) :: square
real(RK) :: simTime = 0.1 ! no. of characteristic time-intervals to simulate
252 character(len=200) :: configFile = "config_file_formatted.in", &
outputFile = "simulation_final_temp_field.dat"
255 call square%init( configFile, simTime ) ! call Initializer
256 call square%run()
258 call square%writeAscii( outputFile )
pause
259 end program solve_heat_diffusion_v1
