C ************************** C * PDE2D 9.2 MAIN PROGRAM * C ************************** C *** 1D PROBLEM SOLVED (GALERKIN METHOD) *** C############################################################################## C Is double precision mode to be used? Double precision is recommended # C on 32-bit computers. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If double precision mode is used, variables and functions assigned +# C + names beginning with a letter in the range A-H or O-Z will be DOUBLE +# C + PRECISION, and you should use double precision constants and FORTRAN +# C + expressions throughout; otherwise such variables and functions will +# C + be of type REAL. In either case, variables and functions assigned +# C + names beginning with I,J,K,L,M or N will be of INTEGER type. +# C + +# C + It is possible to convert a single precision PDE2D program to double +# C + precision after it has been created, using an editor. Just change +# C + all occurrences of "real" to "double precision" +# C + " tdp" to "dtdp" (note leading blank) +# C + Any user-written code or routines must be converted "by hand", of +# C + course. To convert from double to single, reverse the changes. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## implicit double precision (a-h,o-z) parameter (neqnmx= 99) parameter (ndelmx= 20) C############################################################################## C NXGRID = number of X-grid lines # C############################################################################## PARAMETER (NXGRID = 3) C############################################################################## C How many differential equations (NEQN) are there in your problem? # C############################################################################## PARAMETER (NEQN = 2) C DIMENSIONS OF WORK ARRAYS C SET TO 1 FOR AUTOMATIC ALLOCATION PARAMETER (IRWK8Z= 1) PARAMETER (IIWK8Z= 1) PARAMETER (NXP8Z=1001,KDEG8Z=1) C############################################################################## C The solution is saved on a uniform grid of NX+1 points # C XA + I*(XB-XA)/NX # C I=0,...,NX. Enter a value for NX (suggested value = 50). # C # C############################################################################## PARAMETER (NX = 100) PARAMETER (NSAVE = 1) common/parm8z/ pi,D ,RL ,Q dimension xgrid(nxgrid),ixarc(2),xout8z(0:nx),xcross(100),tout8z(0 &:nsave),uout8z(0:nx,2*neqn,0:nsave),uout(0:nx,2,neqn,0:nsave),xres &8z(nxp8z),ures8z(neqn,nxp8z) equivalence (uout,uout8z) allocatable iwrk8z(:),rwrk8z(:) C dimension iwrk8z(iiwk8z),rwrk8z(irwk8z) character*40 title logical plot,symm,fdiff,evcmpx,crankn,fixl,fixr,noupdt,adapt,econ8 &z,ncon8z,restrt,gridid common/dtdp14/ sint8z(20),bint8z(20),slim8z(20),blim8z(20) common/dtdp15/ evlr8z,ev0r,evli8z,ev0i,evcmpx common/dtdp16/ p8z,evr8z(50),evi8z(50) common/dtdp19/ toler(neqnmx),adapt common/dtdp20/ xd0(ndelmx),ndel,idel8z,jdel8z common/dtdp42/ nxa8z,kd8z common/dtdp43/ work8z(nxp8z+3) common/dtdp30/ econ8z,ncon8z common/dtdp46/ eps8z,cgtl8z,npmx8z,itype common/dtdp62/ amin8z(2*neqnmx),amax8z(2*neqnmx) pi = 4.0*atan(1.d0) nxa8z = nxp8z kd8z = kdeg8z C############################################################################## C If you don't want to read the FINE PRINT, default NPROB. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you want to solve several similar problems in the same run, set +# C + NPROB equal to the number of problems you want to solve. Then NPROB +# C + loops through the main program will be done, with IPROB=1,...,NPROB, +# C + and you can make the problem parameters vary with IPROB. NPROB +# C + defaults to 1. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## NPROB = 1 do 78755 iprob=1,nprob C############################################################################## C You may now define global parameters, which may be referenced in any # C of the "FORTRAN expressions" you input throughout the rest of this # C interactive session. You will be prompted alternately for parameter # C names and their values; enter a blank name when you are finished. # C # C Parameter names are valid FORTRAN variable names, starting in # C column 1. Thus each name consists of 1 to 6 alphanumeric characters, # C the first of which must be a letter. If the first letter is in the # C range I-N, the parameter must be an integer. # C # C Parameter values are either FORTRAN constants or FORTRAN expressions # C involving only constants and global parameters defined on earlier # C lines. They may also be functions of the problem number IPROB, if # C you are solving several similar problems in one run (NPROB > 1). Note # C that you are defining global CONSTANTS, not functions; e.g., parameter # C values may not reference any of the independent or dependent variables # C of your problem. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you define other parameters here later, using an editor, you must +# C + add them to COMMON block /PARM8Z/ everywhere this block appears, if +# C + they are to be "global" parameters. +# C + +# C + The variable PI is already included as a global parameter, with an +# C + accurate value 3.14159... +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## D = & 10 RL = & 2.5 Q = & -1 C############################################################################## C A Galerkin finite element method is used, with piecewise polynomial # C basis functions on the subintervals defined by the grid points: # C XGRID(1),XGRID(2),...,XGRID(NXGRID) # C You will first be prompted for NXGRID, the number of X-grid points, # C then for XGRID(1),...,XGRID(NXGRID). Any points defaulted will be # C uniformly spaced between the points you define; the first and last # C points cannot be defaulted. The interval over which the PDE system # C is to be solved is then: # C XGRID(1) < X < XGRID(NXGRID) # C # C############################################################################## call dtdpwx(xgrid,nxgrid,0) C XGRID DEFINED XGRID(1) = & 0 XGRID(NXGRID) = & RL call dtdpwx(xgrid,nxgrid,1) C############################################################################## C Enter the polynomial degree (1,2,3 or 4) desired. A suggested value is # C IDEG = 3. # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + The spatial discretization error is O(h**2), O(h**3), O(h**4) or +# C + O(h**5) when polynomials of degree 1,2,3 or 4, respectively, are +# C + used. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## IDEG = 4 C *******STEADY-STATE PROBLEM itype = 1 t0 = 0.0 dt = 1.0 crankn = .false. noupdt = .false. C Number of Newton iterations NSTEPS = 1 C############################################################################## C If you don't want to read the FINE PRINT, enter 'no'. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Is this problem symmetric? The problem is called symmetric if the +# C + Jacobian matrix +# C + +# C + F1.U F1.Ux ... +# C + A1.U A1.Ux ... +# C + . . +# C + . . +# C + is always symmetric, where F1.U means d(F1)/d(U), and similarly +# C + for the other terms. The Jacobian of GB must also be symmetric, +# C + that is, d(GB_i)/d(U_j) = d(GB_j)/d(U_i), and GB_i must not depend +# C + on Ux,... +# C + +# C + The memory and execution time are decreased if the problem is known +# C + to be symmetric. If you answer 'yes' and your problem is not +# C + symmetric, you will usually get a warning message, but you may get +# C + incorrect answers with no warning, so only answer 'yes' if you are +# C + sure. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## SYMM = .FALSE. FDIFF = .FALSE. C############################################################################## C You may calculate one or more integrals (over the entire region) of # C some functions of the solution and its derivatives. How many integrals # C (NINT), if any, do you want to calculate? # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + In the FORTRAN program created by the preprocessor, the computed +# C + values of the integrals will be returned in the vector SINT8Z. If +# C + several iterations or time steps are done, only the last computed +# C + values are saved in SINT8Z (all values are printed). +# C + +# C + A limiting value, SLIM8Z(I), for the I-th integral can be set +# C + below in the main program. The computations will then stop +# C + gracefully whenever SINT8Z(I) > SLIM8Z(I), for any I=1...NINT. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## NINT = 1 C############################################################################## C You may calculate one or more boundary integrals (over the entire # C boundary) of some functions of the solution and its derivatives. How # C many boundary integrals (NBINT), if any, do you want to calculate? # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + In the FORTRAN program created by the preprocessor, the computed +# C + values of the integrals will be returned in the vector BINT8Z. If +# C + several iterations or time steps are done, only the last computed +# C + values are saved in BINT8Z (all values are printed). +# C + +# C + A limiting value, BLIM8Z(I), for the I-th boundary integral can be +# C + set below in the main program. The computations will then stop +# C + gracefully whenever BINT8Z(I) > BLIM8Z(I), for any I=1...NBINT. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## NBINT = 0 ndel = 0 RESTRT = .FALSE. GRIDID = .TRUE. C############################################################################## C Do you have 'fixed' boundary conditions at the LEFT endpoint, that # C is, are you going to specify values for all unknowns there? # C If you answer 'no', you can specify boundary conditions of the # C more general 'free' type. # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If some, but not all, of the unknowns are to be specified, you +# C + should specify 'free' boundary conditions, since even 'fixed' +# C + type conditions of the form: +# C + Ui = FBi(X,[T]) +# C + can be expressed as 'free' boundary conditions in the form: +# C + Ai*nx = zero(Ui-FBi(X,[T])) +# C + where zero(f) = Big_Number*f is a PDE2D-supplied function. +# C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++# C + Note: while 'fixed' boundary conditions are enforced exactly, 'free' +# C + boundary conditions are not; the greater the overall solution +# C + accuracy, the more closely they are satisfied. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## FIXL = .FALSE. C ixarc(1) = 1 if (fixl) ixarc(1) = -1 C############################################################################## C Do you have 'fixed' boundary conditions at the RIGHT endpoint, that # C is, are you going to specify values for all unknowns there? # C If you answer 'no', you can specify boundary conditions of the # C more general 'free' type. # C############################################################################## FIXR = .FALSE. C ixarc(2) = 2 if (fixr) ixarc(2) = -2 C############################################################################## C The solution is saved on a uniform grid of NX+1 points, covering the # C interval (XA,XB). Enter values for XA,XB. These variables are usually # C defaulted. # C # C The defaults are XA = XGRID(1), XB = XGRID(NXGRID) # C # C############################################################################## C defaults for xa,xb xa = xgrid(1) xb = xgrid(nxgrid) C DEFINE XA,XB HERE: call dtdpx1(nx,xa,xb,hx8z,xout8z,npts8z) C SOLUTION SAVED EVERY NOUT ITERATIONS NOUT = NSTEPS C *******allocate workspace call dtdpzg(nxgrid,ideg,symm,neqn,ii8z,ir8z) if (iiwk8z.gt.1) ii8z = iiwk8z if (irwk8z.gt.1) ir8z = irwk8z allocate (iwrk8z(ii8z),rwrk8z(ir8z)) C *******DRAW GRID POINTS? PLOT = .FALSE. C *******call pde solver call dtdpgx(xgrid, nxgrid, ixarc, restrt, gridid, neqn, ideg, nste &ps, nout, t0, dt, plot, symm, fdiff, itype, nint, nbint, crankn, n &oupdt, xout8z, uout, npts8z, tout8z, nsave, iwrk8z, ii8z, rwrk8z, &ir8z) deallocate (iwrk8z,rwrk8z) C *******read from restart file to array ures8z C call dtdpr1(1,xres8z,nxp8z,ures8z,neqn) C *******write array ures8z back to restart file C call dtdpr1(2,xres8z,nxp8z,ures8z,neqn) C *******call user-written postprocessor call postpr(tout8z,nsave,xout8z,nx,uout,neqn) C *******LINE PLOTS C############################################################################## C Enter a value for IVAR, to select the variable to be plotted or # C printed: # C IVAR = 1 means U (possibly as modified by UPRINT,..) # C 2 A1 # C 3 RM # C 4 A2 # C############################################################################## IVAR = 1 C X IS VARIABLE ics8z = 1 ISET1 = 1 ISET2 = NSAVE ISINC = 1 C alow = amin8z(ivar) ahigh = amax8z(ivar) C############################################################################## C Specify the range (UMIN,UMAX) for the dependent variable axis. UMIN # C and UMAX are often defaulted. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + By default, each plot will be scaled to just fit in the plot area. +# C + For a common scaling, you may want to set UMIN=ALOW, UMAX=AHIGH. +# C + ALOW and AHIGH are the minimum and maximum values over all output +# C + points and over all saved time steps or iterations. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## UMIN = 0.0 UMAX = 0.0 C############################################################################## C Enter a title, WITHOUT quotation marks. A maximum of 40 characters # C are allowed. The default is no title. # C############################################################################## TITLE = ' ' TITLE = 'BEAM PROBLEMS ' call dtdprx(tout8z,nsave,iset1,iset2,isinc) do 78756 is8z=iset1,iset2,isinc call dtdpzp(ics8z,ivar,tout8z,nsave,xout8z,nx,uout8z,neqn,title,um &in,umax,ix8z,is8z) 78756 continue 78755 continue call endgks stop end subroutine pdes8z (yd8z,i8z,j8z,kint8z,x,t) implicit double precision (a-h,o-z) parameter (neqnmx= 99) parameter (ndelmx= 20) C un8z(1,I),un8z(2,I) hold the (rarely used) values C of UI,UIx from the previous iteration or time step common/dtdp48/un8z(2,neqnmx),uu8z(2,neqnmx) common/dtdp11/normx common/dtdp20/xd0(ndelmx),ndel,idel8z,jdel8z double precision normx,delamp(ndelmx,neqnmx) common/parm8z/ pi,D ,RL ,Q U = uu8z(1, 1) Ux = uu8z(2, 1) RM = uu8z(1, 2) RMx = uu8z(2, 2) if (i8z.eq.0) then yd8z = 0.0 C############################################################################## C Enter FORTRAN expressions for the functions whose integrals are to be # C calculated and printed. They may be functions of # C # C X,U,Ux,RM,RMx and (if applicable) T # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you only want to integrate a function over part of the interval, +# C + define that function to be zero on the rest of the interval. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C INTEGRAL DEFINED Utrue = Q/D*(x**4/24 - RL*x**3/6 + RL**2*x**2/4) if (kint8z.eq. 1) yd8z = & abs(U-Utrue) C############################################################################## C Enter FORTRAN expressions for the functions whose "integrals" (sum # C over two boundary points) are to be calculated and printed. They may # C be functions of # C # C X,U,Ux,RM,RMx and (if applicable) T # C # C The unit outward normal, NORMx (=1 at right end point, -1 at left), # C may also be referenced. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you only want to "integrate" a function over one boundary point, +# C + define that function to be zero at the other point. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C BND. INTEGRAL1 DEFINED C if (kint8z.eq.-1) yd8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] else C############################################################################## C Now enter FORTRAN expressions to define the PDE coefficients, which # C may be functions of # C # C X,U,Ux,RM,RMx # C # C and, in some cases, of the parameter T. # C # C Recall that the PDEs have the form # C # C d/dX*A1 = F1 # C d/dX*A2 = F2 # C # C############################################################################## if (j8z.eq.0) then yd8z = 0.0 C F1 DEFINED if (i8z.eq. 1) yd8z = & RM/D C A1 DEFINED if (i8z.eq. 2) yd8z = & Ux C F2 DEFINED if (i8z.eq. 3) yd8z = & Q C A2 DEFINED if (i8z.eq. 4) yd8z = & RMx else endif endif return end function u8z(i8z,x,t0) implicit double precision (a-h,o-z) common/parm8z/ pi,D ,RL ,Q u8z = 0.0 return end function fb8z(i8z,iarc8z,x,t) implicit double precision (a-h,o-z) common/parm8z/ pi,D ,RL ,Q fb8z = 0.0 if (iarc8z.eq.-1) then C FIXL=.FALSE., SO NO FIXED BCs AT LEFT ENDPOINT. C IF FIXL CHANGED TO .TRUE., USE MODEL BELOW TO DEFINE BCs C############################################################################## C Enter FORTRAN expressions to define FB1,FB2 at this endpoint. # C They may be functions of X and (if applicable) T. # C # C Recall that fixed boundary conditions have the form # C # C U = FB1 # C RM = FB2 # C # C############################################################################## C FB1 DEFINED C if (i8z.eq. 1) fb8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] C FB2 DEFINED C if (i8z.eq. 2) fb8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] endif if (iarc8z.eq.-2) then C FIXR=.FALSE., SO NO FIXED BCs AT RIGHT ENDPOINT. C IF FIXR CHANGED TO .TRUE., USE MODEL BELOW TO DEFINE BCs C############################################################################## C Enter FORTRAN expressions to define FB1,FB2 at this endpoint. # C They may be functions of X and (if applicable) T. # C # C Recall that fixed boundary conditions have the form # C # C U = FB1 # C RM = FB2 # C # C############################################################################## C FB1 DEFINED C if (i8z.eq. 1) fb8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] C FB2 DEFINED C if (i8z.eq. 2) fb8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] endif return end subroutine gb8z(gd8z,i8z,j8z,iarc8z,x,t) implicit double precision (a-h,o-z) parameter (neqnmx= 99) C un8z(1,I),un8z(2,I) hold the (rarely used) values C of UI,UIx from the previous iteration or time step. common/dtdp48/un8z(2,neqnmx),uu8z(2,neqnmx) common/dtdp49/bign8z common/parm8z/ pi,D ,RL ,Q zero(f8z) = bign8z*f8z U = uu8z(1, 1) Ux = uu8z(2, 1) RM = uu8z(1, 2) RMx = uu8z(2, 2) if (j8z.eq.0) gd8z = 0.0 if (iarc8z.eq.1) then C############################################################################## C # C Define boundary conditions at the LEFT endpoint X = XGRID(1). # C############################################################################## C############################################################################## C Enter FORTRAN expressions to define GB1,GB2 at this endpoint. They # C may be functions of # C # C X,U,Ux,RM,RMx and (if applicable) T # C # C Recall that free boundary conditions have the form # C # C A1*nx = GB1 (nx=-1 at left endpoint, nx=+1 at right) # C A2*nx = GB2 # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Note that 'fixed' boundary conditions: +# C + Ui = FBi(X,[T]) +# C + can be expressed as 'free' boundary conditions in the form: +# C + GBi = zero(Ui-FBi(X,[T])) +# C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++# C + If Ai=0 (eg, if a PDE is first order), then setting GBi=0 is +# C + equivalent to setting "no" boundary condition (which is sometimes +# C + appropriate for first order PDEs), because the boundary condition +# C + Ai*nx = GBi reduces to 0=0. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C # C############################################################################## if (j8z.eq.0) then C GB1 DEFINED if (i8z.eq. 1) gd8z = & 0 C GB2 DEFINED if (i8z.eq. 2) gd8z = & zero(U) else endif endif if (iarc8z.eq.2) then C############################################################################## C # C Define boundary conditions at the RIGHT endpoint X = XGRID(NXGRID). # C############################################################################## C############################################################################## C Enter FORTRAN expressions to define GB1,GB2 at this endpoint. They # C may be functions of # C # C X,U,Ux,RM,RMx and (if applicable) T # C # C Recall that free boundary conditions have the form # C # C A1*nx = GB1 (nx=-1 at left endpoint, nx=+1 at right) # C A2*nx = GB2 # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Note that 'fixed' boundary conditions: +# C + Ui = FBi(X,[T]) +# C + can be expressed as 'free' boundary conditions in the form: +# C + GBi = zero(Ui-FBi(X,[T])) +# C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++# C + If Ai=0 (eg, if a PDE is first order), then setting GBi=0 is +# C + equivalent to setting "no" boundary condition (which is sometimes +# C + appropriate for first order PDEs), because the boundary condition +# C + Ai*nx = GBi reduces to 0=0. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C # C############################################################################## if (j8z.eq.0) then C GB1 DEFINED if (i8z.eq. 1) gd8z = & zero(RM) C GB2 DEFINED if (i8z.eq. 2) gd8z = & 0 else endif endif return end subroutine pmod8z(x,t,a) implicit double precision (a-h,o-z) parameter (neqnmx= 99) common/dtdp48/un8z(2,neqnmx),uu8z(2,neqnmx) common/dtdp68/upr8z(2,neqnmx),uab8z(2,neqnmx) common/dtdp9/uprint(neqnmx),aprint(neqnmx),bpr8z(neqnmx) common/dtdp14/sint(20),bint(20),slim8z(20),blim8z(20) common/parm8z/ pi,D ,RL ,Q U = uu8z(1, 1) Ux = uu8z(2, 1) a1 = upr8z(2, 1) RM = uu8z(1, 2) RMx = uu8z(2, 2) a2 = upr8z(2, 2) C############################################################################## C If you don't want to read the FINE PRINT, default all of the following # C variables. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Normally, PDE2D saves the values of U,A1,RM,A2 at the output +# C + points. If different variables are to be saved (for later printing +# C + or plotting) the following functions can be used to re-define the +# C + output variables: +# C + define UPRINT(1) to replace U +# C + APRINT(1) A1 +# C + UPRINT(2) RM +# C + APRINT(2) A2 +# C + Each function may be a function of +# C + +# C + X,U,Ux,A1,RM,RMx,A2 and (if applicable) T +# C + +# C + Each may also be a function of the initial triangle number KTRI and +# C + the integral estimates SINT(1),...,BINT(1),... +# C + +# C + The default for each variable is no change, for example, UPRINT(1) +# C + defaults to U. Enter FORTRAN expressions for each of the +# C + following functions (or default). +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C DEFINE UPRINT(*),APRINT(*) HERE: return end C dummy routines function axis8z(i8z,x,y,z,ical8z) implicit double precision (a-h,o-z) axis8z = 0 return end subroutine tran8z(p,q,x,y) implicit double precision (a-h,o-z) return end subroutine xy8z(i8z,iarc8z,s,x,y,s0,sf) implicit double precision (a-h,o-z) return end subroutine dis8z(x,y,ktri,triden,shape) implicit double precision (a-h,o-z) return end subroutine postpr(tout,nsave,xout,nx,uout,neqn) implicit double precision (a-h,o-z) dimension xout(0:nx),tout(0:nsave) dimension uout(0:nx,2,neqn,0:nsave) common/parm8z/ pi,D ,RL ,Q common /dtdp27/ itask,npes,icomm common /dtdp46/ eps8z,cgtl8z,npmx8z,itype data lun,lud/0,47/ if (itask.gt.0) return C UOUT(I,IDER,IEQ,L) = U-sub-IEQ, if IDER=1 C A-sub-IEQ, if IDER=2 C (possibly as modified by UPRINT,..) C at the point XOUT(I) C at time/iteration TOUT(L). C ******* ADD POSTPROCESSING CODE HERE: C IN THE EXAMPLE BELOW, MATLAB PLOTFILES pde2d.m, C pde2d.rdm CREATED (REMOVE C! COMMENTS TO ACTIVATE) C! if (lun.eq.0) then C! lun = 46 C! open (lun,file='pde2d.m') C! open (lud,file='pde2d.rdm') C! endif C! do 78753 l=0,nsave C! if (tout(l).ne.dtdplx(2)) nsave0 = l C!78753 continue C! write (lud,78754) nsave0 C! write (lud,78754) neqn C! write (lud,78754) nx C!78754 format (i8) C! do 78755 i=0,nx C! write (lud,78760) xout(i) C!78755 continue C! do 78759 l=0,nsave0 C! write (lud,78760) tout(l) C! do 78758 ieq=1,neqn C! do 78757 ider=1,2 C! do 78756 i=0,nx C! write (lud,78760) uout(i,ider,ieq,l) C!78756 continue C!78757 continue C!78758 continue C!78759 continue C!78760 format (e16.8) C! write (lun,*) '% Read solution from pde2d.rdm' C! write (lun,*) 'fid = fopen(''pde2d.rdm'');' C! write (lun,*) 'NSAVE = fscanf(fid,''%g'',1);' C! write (lun,*) 'NEQN = fscanf(fid,''%g'',1);' C! write (lun,*) 'NX = fscanf(fid,''%g'',1);' C! if (itype.eq.2) then C! write (lun,*) 'L0 = 0;' C! else C! write (lun,*) 'L0 = 1;' C! endif C! write (lun,*) 'T = zeros(NSAVE+1,1);' C! write (lun,*) 'X = zeros(NX+1,1);' C! write (lun,*) 'U = zeros(NX+1,NSAVE+1,2,NEQN);' C! write (lun,*) 'for i=0:NX' C! write (lun,*) ' X(i+1) = fscanf(fid,''%g'',1);' C! write (lun,*) 'end' C! write (lun,*) 'for l=0:NSAVE' C! write (lun,*) 'T(l+1) = fscanf(fid,''%g'',1);' C! write (lun,*) 'for ieq=1:NEQN' C! write (lun,*) 'for ider=1:2' C! write (lun,*) 'for i=0:NX' C! write (lun,*) C! & ' U(i+1,l+1,ider,ieq) = fscanf(fid,''%g'',1);' C! write (lun,*) 'end' C! write (lun,*) 'end' C! write (lun,*) 'end' C! write (lun,*) 'end' C! write (lun,*) 'xmin = min(X(:));' C! write (lun,*) 'xmax = max(X(:));' C! write (lun,*) '% Plots of each variable' C! write (lun,*) 'for IEQ=1:NEQN' C! write (lun,*) 'IDER = 1;' C! write (lun,*) 'umin = min(min(U(:,L0+1:NSAVE+1,IDER,IEQ)));' C! write (lun,*) 'umax = max(max(U(:,L0+1:NSAVE+1,IDER,IEQ)));' C! write (lun,*) 'for L=L0:NSAVE' C! write (lun,*) ' figure' C! write (lun,*) ' plot(X,U(:,L+1,IDER,IEQ))' C! write (lun,*) ' axis([xmin xmax umin umax])' C! write (lun,*) ' xlabel(''X'')' C! write (lun,*) ' ylabel([''U'',num2str(IEQ)])' C! write (lun,*) ' title(['' T = '',num2str(T(L+1))])' C! write (lun,*) 'end' C! write (lun,*) 'end' return end