function uexa(x,t,kappa,l,io,task) c General information: c Exact solution of the heat equation in an insulated uniform bar of c length l, c u_t = kappa u_xx c with a linear initial temperature gradient at time t=0; c u(x,0) = umax x/l c (with umax a constant read in from file n_linear.in), c and homogeneous Neumann boundary conditions at x=0 and x=l: c u_x(0,t) = u_x(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 'n_linear.in' of c the form: c c c c where umax is the initial maximum temperature in the bar. c For better readability, comment lines may be inserted anywhere c in file n_linear.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 1+2=3 to indicate that the solution is symmetric 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 syml,symr parameter (syml=1,symr=2) 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 cosine 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=syml+symr 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', & '- linear initial temperature distribution with',io,0) call cwrite('uexa', & ' zero temperature at the left end of the bar and',io,0) call cwrite('uexa', & ' maximum temperature at the right end of the bar;',io,0) call cwrite('uexa', & ' insulated ends;',io,0) call cwrite('uexa', & '- n_linear/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','n_linear.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','n_linear.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 t = 0 or kappa = 0, return the linear initial profile: if(t.eq.zero .or. kappa.eq.zero)then uexa=umax*(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 Initialize the temperature to the average: uexa=o2*umax c Sum up to klim terms: do 100 k=0,klim c Evaluate the argument of the cosine, 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 cos(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 -cos(args*x) factor: tk=four*umax/(pi*(two*rk+one))**2*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*cos(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