program bt2cs c General information: c Implicit solution of unsteady heat conduction in a uniform bar. c The equation solved is c u_t = kappa u_xx c for 0 < x < l and 0 < t < tmax c where u is the temperature, l is the length of the bar and tmax is c the desired final time to be computed. c The numerical scheme is the two-step second-order backward-time c scheme. The first step, the first order backward time scheme is c used instead. c Copyright 1996 Leon van Dommelen c Version 1.0 Leon van Dommelen 12/16/96 c Usage information: c Program bt2cs looks for its parameters in a file 'bt2cs.in' of the form: c c c c c c c c where kappa, l and tmax are as above and jmax is the number of mesh c intervals to use along the bar. Value umag should be the typical c magnitude of the solution u to be computed. Program bt2cs needs c a ball-park number in its handling of the boundary conditions. c For better readability, comment lines may be inserted anywhere c in file bt2cs.in. All lines are considered comments, unless the c first nonblank character is a digit, sign, or a decimal point. c During execution, the program will prompt for a numerical time c step to use. It will list some possible choices for this time c step. The program may slightly adjust the chosen time step to c exactly reproduce the desired ending time tmax. c Program output is written to output file 'bt2cs.out'. Note that c the output subroutine out may define additional output files. c The way bt2cs handles the boundary conditions is of note. Unless the c boundary values of u are constant, program bt2cs will attempt to c satisfy the boundary conditions by a process of trial and error. c More precisely, program bt2cs will perturb the left and right c boundary values by selected amounts, and then compute fractions of c these perturbations that satisfy the boundary conditions as c returned by subroutines bc1 and bc2. c External variables and info for compiling or changing program bt2cs: c The next statement can be removed if your compiler does not support it: implicit none c Program bt2cs requires a separate subroutine: c initx(l,jdim,jmax,ndim,io,x,t,u) c to determine the initial mesh point locations x(j). c Subroutine initx can be found in file: c ../../init/init_x.f external initx c Program bt2cs requires a separate subroutine: c initu(kappa,l,jdim,jmax,ndim,n,x,t,io,task,u) c to set the initial values for the temperature u(j,0). c For an example subroutine initu see: c ../../init/exact/init_u.f c For a template for writing new subroutines initu see: c ../../init/template/init_u.f external initu c Program bt2cs requires a separate subroutine: c bc1(kappa,l,jdim,jmax,ndim,n,x,t,io,task,u) c to set the left hand boundary condition. c For an example subroutine bc1 see: c ../../bc_1/diri_exa/bc_1.f c For a template for writing new subroutines bc1 see: c ../../bc_1/template/bc_1.f external bc1 c Program bt2cs requires a separate subroutine: c bc2(kappa,l,jdim,jmax,ndim,n,x,t,io,task,u) c to set the right hand boundary condition. c For an example subroutine bc2 see: c ../../bc_2/diri_exa/bc_2.f c For a template for writing new subroutines bc2 see: c ../../bc_2/template/bc_2.f external bc2 c Program bt2cs requires a separate subroutine: c out(kappa,l,jdim,jmax,ndim,n,x,t,u,io,task) c to handle output of the results. c For an example subroutine out see: c ../../out/coef_exa/out.f c For a template for writing new subroutines out see: c ../../out/template/out.f external out c Program bt2cs requires a separate subroutine: c tridia(ndim,n,a,b,c,d,w1,w2,w3,w4,w5,u,ier,err,umax) c to solve a tridiagonal system of equations of the form c b(0) u(0) + c(0) u(1) = d(0) first equation; c a(j) u(j-1) + b(j) u(j) + c(j) u(j+1) = d(j) for j = 0, 1, ..., n-1; c a(n) u(n-1) + b(n) u(n) = d(n) final equation; c for u(0), u(1), ... u(n). Arguments w1 through w5 are work space c and ier is an error code which should be less than two after tridia. c Subroutine tridia stores its estimated error in err and the maximum c absolute value of the u(i) in umax. c Subroutine tridia can be found in: c ../../../lib/tridia.f external tridia c Program bt2cs requires a separate subroutine: c satis2(ndim,n,u0,r10,r20,u1,r11,r21,u2,r12,r22, C erru,errr,u,ier,err) c to combine three vectors u0, u1, and u2 so that two equations c become satisfied. Here r10, r11, and r12 are the residuals in the c first equation for u0, u1, and u2 respectively; r20, r21, and r22 c are the residuals in the second equation; erru is the estimated c maximum absolute error in the vectors ui; errr is the estimated c maximum absolute error in the residuals, u is the combined c solution; ier is an error which should ne nozero after satis2; and c err is the estimated maximum absolute error in u. c Subroutine satis2 can be found in: c ../../../lib/satis2.f external satis2 c Program bt2cs uses the following utility routines from c ../../../lib/util.f: c exit0(module) exits the program after successful completion. external exit0 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 fcreat(module,filnam,io) creates the new output file "filnam", c opening it on I/O unit "io". It selects "io" if initialized to zero. external fcreat 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 iread(module,text,ioin,io,show) reads integer "iread" 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. integer iread external iread c iwrite(module,text,val,io,show) writes integer "val" and its c description "text" to I/O unit "io" if nonzero, and to the screen c if "show" is nonzero. external iwrite c rask(module,text,io) gets the number "rask" from the user, c using the prompt "text", and writes it to I/O unit "io" if nonzero. double precision rask external rask 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 Program variables: c The heat conduction coefficient: double precision kappa c The desired maximum time to compute to: double precision tmax c The size of the time step: double precision dt c The time step number n: integer n c The maximum value of n, corresponding to time tmax: integer nmax c Program bt2cs declares only storage for three time levels. This c is achieved by storing using index nn=mod(n,3) instead of index n c itself. (Note that mod(n,3) equals 0, 1, or 2 depending on the c value of n.) c The value of nn: integer nn c The value of nn corresponding to n+1: integer nnp1 c The value of nn corresponding to n-1: integer nnm1 c The array of mesh times, with only storage for 3 time levels (see c the comments above): double precision t(0:2) c The length of the bar: double precision l c The size of the spatial step: double precision dx c The mesh point index: integer j c The maximum mesh point index: integer jmax c The dimensioned maximum mesh point index: integer jdim parameter (jdim=128) c The x-values of the mesh points: double precision x(0:jdim) c The values x(j) are set by subroutine initx, which will take c the implementation of the boundary conditions into account. The c value of x(0) will be approximately zero, depending on the left c boundary condition, and the value of x(jmax) will be approximately c l, depending on the right boundary condition. c The temperature values of the mesh points, with only storage for c 3 time levels (see the comments above at the declaration of nn): double precision u(0:jdim,0:2) c The Fortran I/O unit number for the output file: integer io c The Fortran I/O unit number for the input file: integer ioin c The coefficients in the discretization of u_t double precision utn,utnp1,utnm1 c The properties of boundary condition subroutines bc1 and bc2: integer ibc1,ibc2 c The typical magnitude of u: double precision umag c The typical amount to perturb the boundary conditions while c trying to satisfy them: double precision pertb c The amount to perturb the left boundary condition: double precision pertb1 c The amount to perturb the right boundary condition: double precision pertb2 c Storage for the solution of the tridiagonal equations: c The coefficients of the tridiagonal equations: double precision a(0:jdim),b(0:jdim),c(0:jdim),d(0:jdim) c The solution with assumed boundary conditions: double precision u0(0:jdim) c The residuals in the two boundary conditions for solution u0: double precision r1u0,r2u0 c The solution with a perturbed left boundary: double precision u1(0:jdim) c The residuals in the two boundary conditions for solution u1: double precision r1u1,r2u1 c The solution with a perturbed right boundary: double precision u2(0:jdim) c The residuals in the two boundary conditions for solution u2: double precision r1u2,r2u2 c The combination of u0, u1, and u2 that satisfies the boundary conditions: double precision uu(0:jdim) c Work space for use by subroutine tridia double precision & w1(0:jdim),w2(0:jdim),w3(0:jdim),w4(0:jdim),w5(0:jdim) c Error code returned by subroutine tridia integer ier c The maximum error in the u(i,n) due to subroutine tridia alone: double precision trierr c The maximum value of |u(i,n)| computed by tridia double precision umax c The estimated error in the residuals: double precision reserr c The estimated error due to solving the implicit equations alone: double precision solerr c An error: double precision err c The machine epsilon (the relative error in real numbers): double precision eps save eps c Various constants, defined to make conversion of precision simpler: double precision zero,one,two,three,four parameter (zero=0.d0,one=1.d0,two=2.d0,three=3.d0,four=4.d0) double precision o2 parameter (o2=.5d0) c Executable statements: c Find the machine epsilon as the absolute error in the value 1: eps=two 10 eps=eps/two if(one+eps.ne.one)goto 10 c Open the output file bt2cs.out on I/O unit io using fcreat of util.f: io=0 call fcreat('bt2cs','bt2cs.out',io) c Open the input file bt2cs.in on I/O unit ioin using fopen of util.f: ioin=0 call fopen('bt2cs','bt2cs.in',ioin) c Read the input file: c Read the heat conduction coefficient: kappa=rread('bt2cs', &'The heat conduction coefficient kappa',ioin,io,1) c Allow negative kappa (the backward heat equation): if(kappa.lt.zero)call warn('bt2cs',io, &'Negative heat conduction coefficient.',' ',' ') c Read the length of the bar: l=rread('bt2cs','The length of the bar l',ioin,io,1) c Enforce a positive bar length: if(l.le.zero)call fatal('bt2cs', &'The length of the bar must be positive.',' ',' ') c Final time to compute: tmax=rread('bt2cs','The maximum time to compute',ioin,io,1) c Do not allow backward marching: if(tmax.lt.zero)call fatal('bt2cs', &'The maximum time cannot be negative.',' ',' ') c Allow tmax to be zero to compute the initial plane only: if(tmax.eq.zero)call warn('bt2cs',io, &'Only initial data will be computed.',' ',' ') c Read the number of mesh intervals along the bar to use: jmax=iread('bt2cs','The number of mesh intervals',ioin,io,1) c Check whether we have enough storage available: if(jmax.gt.jdim)then call iwrite('bt2cs',' bt2cs: '// & ' Requested number of mesh intervals',jmax,io,1) call iwrite('bt2cs',' '// & 'Dimensioned number of mesh intervals',jdim,io,1) call fatal('bt2cs','The requested number of mesh intervals '// & 'exceeds the program limit.',' ',' ') endif c Less than one interval is not supported: if(jmax.lt.1)call fatal('bt2cs', &'The number of mesh intervals should be at least one.',' ',' ') c Typical magnitude of the solution: umag=rread('bt2cs','Typical magnitude of the solution',ioin,io,1) c The magnitude must be positive: if(umag.le.zero)call fatal('bt2cs', &'The magnitude must be positive.',' ',' ') c Done with input call fclose('bt2cs','bt2cs.in',ioin) c Initialize the program variables: c Let subroutine initx create the mesh point positions x(j): call initx(l,jdim,jmax,2,io,x,t,u) c Compute the mesh size: dx=x(1)-x(0) call rwrite('bt2cs','Mesh size dx',dx,io,0) c Let rask read the time step from the keyboard: if(kappa.ne.zero)then call rwrite('bt2cs','dx^2/kappa',dx**2/kappa,io,1) call rwrite('bt2cs','dx l/kappa',dx*l/kappa,io,1) endif 400 dt=one if(tmax.gt.zero)dt=rask('bt2cs','the time step (0 exits)',io) if(dt.eq.zero)call fatal('bt2cs', &'Execution aborted by the user.',' ',' ') if(dt.lt.zero)then call warn('bt2cs',io, & 'The time step must be positive.', & 'The entered negative value has been refused.',' ') goto 400 endif c Compute the rounded number of time steps: nmax=tmax/dt+o2 if(tmax.le.zero)nmax=0 if(tmax.gt.zero)nmax=max(nmax,1) call iwrite('bt2cs','The total number of time steps',nmax,io,1) c Adjust the time step slightly to give an integer amount of steps: if(nmax.gt.0)then dt=tmax/nmax call rwrite('bt2cs','The actual time step to be used',dt,io,1) endif c Output the amount of work to be done: call iwrite('bt2cs', &'The total number of points to be computed',jmax*nmax,io,1) c Start at time zero: t(0)=zero c Make sure that the various subroutines do not mistakenly think the c other levels have meaningful information for earlier times: t(1)=t(0)+dt t(2)=t(1)+dt c Allow subroutine init to set the initial data u(j,0): call initu(kappa,l,jdim,jmax,2,0,x,t,io,0,u) c Allow subroutine bc1 to set the boundary condition at j=0: call bc1(kappa,l,jdim,jmax,2,0,x,t,io,0,u) c Allow subroutine bc2 to set the boundary condition at j=jmax: call bc2(kappa,l,jdim,jmax,2,0,x,t,io,0,u) c Perform output: call out(kappa,l,jdim,jmax,2,0,x,t,u,io,0) c Initialize the error in the solution procedure: solerr=zero c Select how much to perturb the boundary conditions while trying to c satisfy them. Normally we perturb by the typical magnitude of u, c but if the values returned by bc1 or bc2 vary nonlinearly with c changes in the surrounding u-values, we must use small c perturbations to be able to linearize their behavior: c Normal perturbation: pertb=umag c Question bc1 and bc2 about their properties: ibc1=-1 call bc1(kappa,l,jdim,jmax,2,0,x,t,io,ibc1,u) ibc2=-1 call bc2(kappa,l,jdim,jmax,2,0,x,t,io,ibc2,u) c If either ibc1 or ibc2 is odd, the boundary condition is nonlinear: if(mod(ibc1,2).ne.0 .or. mod(ibc2,2).ne.0)pertb=umag*sqrt(eps) c March the solution in time to the final time tmax: c For n= 0 to nmax-1, compute the time plane n+1: do 5000 n=0,nmax-1 c Use the Fortran mod (remainder) function to reduce n-1, n, and c n+1 to the only two available storage levels in array u: nnm1=mod(n-1,3) if(n.eq.0)nnm1=0 nn=mod(n,3) nnp1=mod(n+1,3) c Find the coefficients of the discretization for the time derivative: if(n.eq.0)then c Use a first-order two-point discretization at n = 0: utnp1=one/dt utn=-one/dt utnm1=zero else c Use a second-order three-point discretization: utnp1=three/(two*dt) utn=-four/(two*dt) utnm1=one/(two*dt) endif c Find the time at step n+1: t(nnp1)=(n+1+zero)/(nmax+zero)*tmax c Written this way, the last time will be exactly tmax despite c any round-off errors. So if the last requested output time c equals tmax, it will not be missed because of round-off. c Compute the coefficients of the tridiagonal equations for c the interior points, as given by the Crank-Nicholson scheme: do 2200 j=1,jmax-1 c Subdiagonal coefficient (= coefficient of u(j-1,n+1)): a(j)=-kappa/dx**2 c Diagonal coefficient (= coefficient of u(j,n+1)): b(j)=utnp1+two*kappa/dx**2 c Superdiagonal coefficient (= coefficient of u(j+1,n+1)): c(j)=a(j) c Right hand side (u(j,n) and u(j,n-1) terms): d(j)=-u(j,nn)*utn-u(j,nnm1)*utnm1 2200 continue c We do not know the equations for the boundary points j = 0 and c j = jmax used by the boundary condition subroutines bc1 and c bc2. For now, let us simply assume that the boundary values c will be unchanged at the next time: c The assumed left boundary condition u(0,n+1) = u(0,n): c Diagonal coefficient: b(0)=one c Superdiagonal coefficient; c(0)=zero c Right hand side: d(0)=u(0,nn) c The assumed right boundary condition u(jmax,n+1) = u(jmax,n): c Subdiagonal coefficient: a(jmax)=zero c Diagonal coefficient: b(jmax)=one c Right hand side: d(jmax)=u(jmax,nn) c Let subroutine tridia solve this system: call tridia(jdim,jmax,a,b,c,d,w1,w2,w3,w4,w5,u0,ier,err,umax) if(ier.gt.1)call fatal('bt2cs', & 'Fatal error while solving the tridiagonal system for u0.', & ' ',' ') trierr=err c Copy the solution found by tridia over into the array u: do 3000 j=0,jmax u(j,nnp1)=u0(j) 3000 continue c Let subroutine bc1 return its u(0,n+1) for this solution: call bc1(kappa,l,jdim,jmax,2,nnp1,x,t,io,0,u) c Compute the residual in the left boundary condition: r1u0=u(0,nnp1)-d(0) c Let subroutine bc2 return its u(jmax,n+1) for this solution: call bc2(kappa,l,jdim,jmax,2,nnp1,x,t,io,0,u) c Compute the residual in the right boundary condition: r2u0=u(jmax,nnp1)-d(jmax) c All done if the assumed boundary conditions are correct: if(r1u0.eq.zero .and. r2u0.eq.zero)goto 4900 c Tentatively perturb the left boundary condition: c Select a perturbation amount: pertb1=max(pertb,abs(r1u0)) if(r1u0.lt.zero)pertb1=-pertb1 c Now perturb the right hand side of the first equation: d(0)=u(0,nn)+pertb1 c Let subroutine tridia solve the perturbed system: call tridia(jdim,jmax,a,b,c,d,w1,w2,w3,w4,w5,u1,ier,err,umax) if(ier.gt.1)call fatal('bt2cs', & 'Fatal error while solving the tridiagonal system for u1.', & ' ',' ') trierr=max(trierr,err) c Copy the perturbed solution found by tridia over into the array u: do 3100 j=0,jmax u(j,nnp1)=u1(j) 3100 continue c Let subroutine bc1 return its u(0,n+1) for this perturbed solution: call bc1(kappa,l,jdim,jmax,2,nnp1,x,t,io,0,u) c Compute the residual in the left boundary condition: r1u1=u(0,nnp1)-d(0) c Let subroutine bc2 return its u(jmax,n+1) for this solution: call bc2(kappa,l,jdim,jmax,2,nnp1,x,t,io,0,u) c Compute the residual in the right boundary condition: r2u1=u(jmax,nnp1)-d(jmax) c Restore the first equation: d(0)=u(0,nn) c Tentatively perturb the right boundary condition: c Select a perturbation amount: pertb2=max(pertb,abs(r2u0)) if(r2u0.lt.zero)pertb2=-pertb2 c Now perturb the right hand side of the last equation: d(jmax)=u(jmax,nn)+pertb2 c Let subroutine tridia solve the perturbed system: call tridia(jdim,jmax,a,b,c,d,w1,w2,w3,w4,w5,u2,ier,err,umax) if(ier.gt.1)call fatal('bt2cs', & 'Fatal error while solving the tridiagonal system for u2.', & ' ',' ') trierr=max(trierr,err) c Copy the solution found by tridia over into the array u: do 3200 j=0,jmax u(j,nnp1)=u2(j) 3200 continue c Let subroutine bc1 return its u(0,n+1) for this solution: call bc1(kappa,l,jdim,jmax,2,nnp1,x,t,io,0,u) c Compute the residual in the left boundary condition: r1u2=u(0,nnp1)-d(0) c Let subroutine bc2 return its u(jmax,n+1) for this solution: call bc2(kappa,l,jdim,jmax,2,nnp1,x,t,io,0,u) c Compute the residual in the right boundary condition: r2u2=u(jmax,nnp1)-d(jmax) c Assume that the error in the resduals is of the order of c the error in the u values: reserr=trierr c Let subroutine satis2 figure out the combination of u0, u1, and c u2 that satisfies both boundary conditions: call satis2(jdim,jmax,u0,r1u0,r2u0,u1,r1u1,r2u1,u2,r1u2,r2u2, & trierr,reserr,uu,ier,err) if(ier.gt.1)call fatal('bt2cs', & 'Unable to satisfy both boundary conditions.', & ' ',' ') solerr=max(solerr,err) c Store this solution into the array u: do 3300 j=0,jmax u(j,nnp1)=uu(j) 3300 continue c Allow subroutine out to perform any required output: 4900 call out(kappa,l,jdim,jmax,2,nnp1,x,t,u,io,0) 5000 continue c Exit: c Write out the estimated error made in solving the equations: call rwrite('bt2cs', &'The estimated error made in solving the equations alone', & solerr,io,1) c Close the output file: call fclose('bt2cs','bt2cs.out',io) c Exit: call exit0('bt2cs') c For picky compilers: stop end