function uexa(x,t,kappa,l,io,task) c General information: c Exact solution of the heat equation in a uniform bar of length l, c u_t = kappa u_xx c with a parabolic initial temperature profile at time t=0; c u(x,0) = 4 umax (1 - x/l) x/l for 0 <= x <= l c (with umax a constant read in from file d_para.in), c and homogeneous Dirichlet boundary conditions at x=0 and x=l: c u(0,t) = u(l,t) = 0 c Copyright 1996 Leon van Dommelen c Version 1.0 Leon van Dommelen 12/16/96 c Usage information: c Function uexa looks for its parameters in a file 'd_para.in' of c the form: c c c c where umax is the initial temperature in the middle of the bar. c For better readability, comment lines may be inserted anywhere c in file d_para.in. All lines are considered comments, unless the c first nonblank character is a digit, sign, or a decimal point. c Function uexa can be called in three different ways depending on c the chosen value of parameter task: c 1. task = 0: c This is the standard call. In this case, uexa computes and c returns the exact temperature at position (x,t). c During such a call, all input parameters must have the c appropriate values. c 2. task = -1: c This is the inquiry call. Function uexa changes the value of c task to indicate its properties. To be precise, its sets task to c 8+4=12 to indicate that the solution is antisymmetric around both c boundaries. c 3. Any other value: c If task is positive, subroutine uexa merely initializes itself c and exits, otherwise it does nothing at all. c Arguments: c Avoid typos: implicit none c Input: position x, (0 <= x <= l), at which the temperature is needed: double precision x c Input: time t, (0 <= t), at which the temperature is needed: double precision t c Input: conduction constant (0 <= kappa): double precision kappa c Input: length of the bar (0 < l): double precision l c Input: I/O unit of an already open output file or zero: integer io c Function uexa only writes to this unit if io is positive. c Input/Output: task to perform: integer task c See the usage information above for more information on task. c Output: The temperature at location x and time t (if task = 0): double precision uexa c If task is nonzero, the returned value uexa is always zero. c External variables and info for compiling or changing function uexa: c The following utility routines from ../../../lib/util.f were used: c cwrite(module,text,io,show) writes line "text" to I/O unit c "io" if nonzero, and to the screen if "show" is nonzero. external cwrite c fatal(module,text1,text2,text3) kills the program after a c fatal error, printing the lines "text1", "text2" and "text3". external fatal c fclose(module,filnam,io) closes the I/O unit "io" on which c file "filnam" is open. external fclose c fopen(module,filnam,io) opens the existing file "filnam" on c I/O unit "io". It selects "io" if initialized to zero. external fopen c rread(module,text,ioin,io,show) reads number "rread" from I/O c unit "ioin". It writes it and its description "text" to I/O unit c "io" if nonzero, and to the screen if "show" is nonzero. double precision rread external rread c rwrite(module,text,val,io,show) writes number "val" and its c description "text" to I/O unit "io" if nonzero, and to the screen c if "show" is nonzero. external rwrite c warn(module,io,text1,text2,text3) writes the warning lines c "text1", "text2" and "text3" to the screen, and to I/O unit "io" c if nonzero. external warn c Local variables: c An integer to keep track of whether uexa has been initialized: integer init c Integers that are returned during the inquiry call (see above): integer asyml,asymr parameter (asyml=4,asymr=8) c Maximum initial temperature: double precision umax save umax c The exact solution can be written as a Fourier series using the c standard separation-of-variable procedure. Function uexa simply c sums that series. Used variables: c Index of the Fourier series: integer k c The same index as a floating point number (used since conversions c from integer to floating point may be very slow): double precision rk c Individual term in the Fourier series: double precision tk c Argument of the sine in the individual term: double precision args c Argument of the exponential factor in the individual term: double precision arge c A value for x above which exp(x) does not underflow (used since c underflow error messages can slow down execution greatly): double precision expufl parameter (expufl=-70.d0) c Maximum number of terms ever to sum in the Fourier series: integer klim parameter (klim=10000000) c Regardless of klim, the sum will terminate when converged. c Function uexa will crash if klim is exceeded. c Fortran I/O unit number for the input file: integer ioin c Constants defined to make changing precision easier: double precision zero,one,two,four,o2 parameter (zero=0.d0,one=1.d0,two=2.d0,four=4.d0,o2=0.5d0) double precision pi save pi c Data statements: data init/0/ c Executable statements: c Default value of the temperature: uexa=zero c For task = -1, return the properties of the solution: if(task.eq.-1)then c Both boundaries are anti-symmetric: task=asyml+asymr return endif c Ignore other negative tasks: if(task.lt.0)return c Initialization during the first time that out is called: if(init.eq.0)then c Show which function is being initialized: call cwrite('uexa', & 'Initialization of exact solution function uexa:',io,0) call cwrite('uexa', & '- parabolic initial temperature distribution with',io,0) call cwrite('uexa', & ' maximum temperature in the center of the bar;',io,0) call cwrite('uexa', & '- zero temperature at the ends of the bar.',io,0) call cwrite('uexa', & '- d_para/u_exa.f Copyright 1996 Leon van Dommelen.',io,0) c Value of pi: pi=four*atan(one) c Open the input file: ioin=0 call fopen('uexa','d_para.in',ioin) c Read in the maximum temperature to use: umax=rread('uexa','The maximum initial temperature',ioin,io,1) c Accept a trivial solution: if(umax.eq.zero)call warn('uexa',io, & 'Zero initial temperature accepted.',' ',' ') c Close the input file: call fclose('uexa','d_para.in',ioin) c All done with the initialization: init=1 endif c Further ignore nonzero tasks: if(task.ne.0)return c Check the arguments: c The Fourier series will not converge for negative kappa: if(kappa.lt.zero)then call rwrite('uexa', & ' uexa: Heat conduction coefficient kappa',kappa,io,1) call fatal('uexa', & 'The heat conduction coefficient must be positive.',' ',' ') endif c Disallow negative length to simplify the code: if(l.le.zero)then call rwrite('uexa',' uexa: Bar length l',l,io,1) call fatal('uexa','The bar length must be positive.',' ',' ') endif c The Fourier series does not converge for negative times: if(t.lt.zero)then call rwrite('uexa',' uexa: Time t',t,io,1) call fatal('uexa','Negative times are not allowed.',' ',' ') endif c Disallow positions outside the bar to simplify the code: if(x.lt.zero .or. x.gt.l)then call rwrite('uexa',' uexa: Bar length l',l,io,1) call rwrite('uexa',' uexa: Position x',x,io,1) call fatal('uexa', & 'Position x should be in the range from 0 to l.',' ',' ') endif c Compute the temperature: c For umax = 0, return uexa as zero: if(umax.eq.zero)goto 900 c For x = 0 or x = l, return the homogenuous boundary condition: if(x.eq.zero .or. x.eq.l)goto 900 c For t = 0 or kappa = 0, return the parabolic initial profile: if(t.eq.zero .or. kappa.eq.zero)then uexa=four*umax*(one-x/l)*x/l goto 900 endif c Else find the exact solution from summing the Fourier series: c Initialize the floating point value of the Fourier series index k: rk=zero c Sum up to klim terms: do 100 k=0,klim c Evaluate the argument of the sine, except for the factor x: args=(two*rk+one)*pi/l c Evaluate the argument of the exponential: arge=-args**2*kappa*t c Stop summing when underflow of the exponential is imminent: if(arge.lt.expufl)goto 900 c The sin(args*x) term can make convergence hard to judge since c it can oscillate. So we will judge convergence based on its c amplitude alone, in the following way: c Evaluate the k-th term except for the sin(args*x) factor: tk=two*four**2*umax/(pi*(two*rk+one))**3*exp(arge) c Terminate when tk is too small to change the current value: if(tk+uexa.eq.uexa)goto 900 c finish the k-th term: tk=tk*sin(args*x) c Add the term to the sum: uexa=uexa+tk c Increment the real value of k and go do the next term: rk=rk+one 100 continue c We failed to get convergence since we did not jump out of the loop: call rwrite('uexa', &' uexa: Last term in the Fourier series',tk,io,1) call rwrite('uexa', &' Current value for the sum',uexa,io,1) call fatal('uexa','The Fourier series did not converge.',' ',' ') c Exit: c Jump here when done: 900 return end