C ************************** C * PDE2D 9.5 MAIN PROGRAM * C ************************** C *** 0D PROBLEM SOLVED *** C############################################################################## C Is double precision mode to be used? Double precision is recommended. # 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) C############################################################################## C How many differential equations (NEQN) are there in your problem? # C############################################################################## PARAMETER (NEQN = 5) parameter (nxgrid = 1) C DIMENSIONS OF WORK ARRAYS C SET TO 1 FOR AUTOMATIC ALLOCATION PARAMETER (IRWK8Z= 1) PARAMETER (IIWK8Z= 1) C############################################################################## C The solution will be saved (for possible postprocessing) at the NSAVE+1 # C time points # C T0 + K*(TF-T0)/NSAVE # C K=0,...,NSAVE. Enter a value for NSAVE. # C # C If a user-specified constant time step is used, NSTEPS must be an # C integer multiple of NSAVE. # C############################################################################## PARAMETER (NSAVE = 200) common/parm8z/ pi,RL ,G ,Rm dimension xgrid(nxgrid),xout8z(0:0),tout8z(0:nsave),uout(2,neqn,0: &nsave) allocatable iwrk8z(:),rwrk8z(:) C dimension iwrk8z(iiwk8z),rwrk8z(irwk8z) character*40 title logical linear,crankn,noupdt,nodist,fillin,evcmpx,adapt,plot,lsqfi &t,fdiff,econ8z,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/dtdp30/ econ8z,ncon8z common/dtdp46/ eps8z,cgtl8z,npmx8z,itype,near8z pi = 4.0*atan(1.d0) 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 PDE2D solves the time-dependent system (note: U,F,U0 may be vectors, # C C,RHO may be matrices): # C # C C(T,U)*d(U)/dT = F(T,U) # C # C or the algebraic system: # C # C F(U) = 0 # C # C or the linear and homogeneous eigenvalue system: # C # C F(U) = lambda*RHO*U # C # C For time-dependent problems there are also initial conditions: # C # C U = U0 at T=T0 # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + A system of NEQN complex equations must be written as a system +# C + 2*NEQN real equations, by separating the equations into their real +# C + and imaginary parts. However, note that the complex arithmetic +# C + abilities of FORTRAN can be used to simplify this separation. For +# C + example, the complex nonlinear algebraic equation: +# C + U**10 - 1.0 = 0, where U = UR + UI*I +# C + would be difficult to split up analytically, but using FORTRAN +# C + expressions it is easy: +# C + F1 = REAL(CMPLX(UR,UI)**10-1.0) +# C + F2 = AIMAG(CMPLX(UR,UI)**10-1.0) +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# 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############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- RL = & 1.0 G = & 32.0 Rm = & 4.0 C xgrid(1) = 0.0 C *******TIME-DEPENDENT PROBLEM itype = 2 C############################################################################## C Enter the initial time value (T0) and the final time value (TF), for # C this time-dependent problem. T0 defaults to 0. # C # C TF is not required to be greater than T0. # C############################################################################## T0 = 0.0 C-----------------------------------------> INPUT FROM GUI <------------------- T0 = & 0 C-----------------------------------------> INPUT FROM GUI <------------------- TF = & 5 C############################################################################## C Is this a linear problem? ("linear" means all differential equations # C and all boundary conditions are linear). If you aren't sure, it is # C safer to answer "no". # C############################################################################## LINEAR = .FALSE. C############################################################################## C Do you want the time step to be chosen adaptively? If you answer # C 'yes', you will then be prompted to enter a value for TOLER(1), the # C local relative time discretization error tolerance. The default is # C TOLER(1)=0.01. If you answer 'no', a user-specified constant time step # C will be used. We suggest that you answer 'yes' and default TOLER(1) # C (although for certain linear problems, a constant time step may be much # C more efficient). # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If a negative value is specified for TOLER(1), then ABS(TOLER(1)) is +# C + taken to be the "absolute" error tolerance. If a system of PDEs is +# C + solved, by default the error tolerance specified in TOLER(1) applies +# C + to all variables, but the error tolerance for the J-th variable can +# C + be set individually by specifying a value for TOLER(J) using an +# C + editor, after the end of the interactive session. +# C + +# C + Each time step, two steps of size dt/2 are taken, and that solution +# C + is compared with the result when one step of size dt is taken. If +# C + the maximum difference between the two answers is less than the +# C + tolerance (for each variable), the time step dt is accepted (and the +# C + next step dt is doubled, if the agreement is "too" good); otherwise +# C + dt is halved and the process is repeated. Note that forcing the +# C + local (one-step) error to be less than the tolerance does not +# C + guarantee that the global (cumulative) error is less than that value.+# C + However, as the tolerance is decreased, the global error should +# C + decrease correspondingly. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- ADAPT = .FALSE. TOLER(1) = 0.01 NOUPDT = .FALSE. C############################################################################## C The time stepsize will be chosen adaptively, between an upper limit # C of DTMAX = (TF-T0)/NSTEPS and a lower limit of 0.0001*DTMAX. Enter # C a value for NSTEPS (the minimum number of steps). # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you later turn off adaptive time step control, the time stepsize +# C + will be constant, DT = (TF-T0)/NSTEPS. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C-----------------------------------------> INPUT FROM GUI <------------------- NSTEPS = & 1000 dt = (tf-t0)/max(nsteps,1) C############################################################################## C If you don't want to read the FINE PRINT, enter 'yes'. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Is the Crank-Nicolson scheme to be used to discretize time? If you +# C + answer 'no', a backward Euler scheme will be used. Do not use the +# C + Crank Nicolson method if the left hand side of any ODE is zero, for +# C + example, if a differential/algebraic system is solved. +# C + +# C + The Crank-Nicolson scheme is second order, and the backward Euler +# C + method is first order. However, if adaptive time step control is +# C + chosen, an extrapolation is done between the 1-step and 2-step +# C + answers which makes the Euler method second order, and the Crank- +# C + Nicolson method third order (fourth order for linear problems). +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## CRANKN = .FALSE. C############################################################################## C If you don't want to read the FINE PRINT, enter 'yes' (strongly # C recommended). # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + The partial derivatives of some of the PDE and boundary condition +# C + coefficients are required by PDE2D. These may be calculated +# C + automatically using a finite difference approximation, or supplied +# C + by the user. Do you want them to be calculated automatically? +# C + +# C + If you answer 'yes', you will not be asked to supply the derivatives,+# C + but there is a small risk that the inaccuracies introduced by the +# C + finite difference approximation may cause the Newton iteration +# C + to converge more slowly or to diverge, especially if low precision +# C + is used. This risk is very low, however, and since answering 'no' +# C + means you may have to compute many partial derivatives, it is +# C + recommended you answer 'yes' unless you have some reason to believe +# C + there is a problem with the finite difference approximations. +# C + +# C + If you supply analytic partial derivatives, PDE2D will do some spot +# C + checking and can usually issue a warning if any are supplied +# C + incorrectly. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## FDIFF = .TRUE. 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 = 0 nbint = 0 lsqfit = .false. C############################################################################## C If you don't want to read the FINE PRINT, enter 'no'. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Do you want to read the initial conditions from the restart file, +# C + if it exists (and use the conditions supplied above if it does not +# C + exist)? +# C + +# C + If so, PDE2D will dump the final solution at the end of each run +# C + into a restart file "pde2d.res". Thus the usual procedure for +# C + using this dump/restart option is to make sure there is no restart +# C + file in your directory left over from a previous job, then the +# C + first time you run this job, the initial conditions supplied above +# C + will be used, but on the second and subsequent runs the restart file +# C + from the previous run will be used to define the initial conditions. +# C + +# C + You can do all the "runs" in one program, by setting NPROB > 1. +# C + Each pass through the DO loop, T0,TF,NSTEPS and possibly other +# C + parameters may be varied, by making them functions of IPROB. +# C + +# C + If the 2D or 3D collocation method is used, the coordinate +# C + transformation should not change between dump and restart. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## RESTRT = .FALSE. gridid = .true. iperdc = 0 npts8z = 1 xout8z(0) = 0.0 call dtdp1q(nxgrid,neqn,ii8z,ir8z) if (iiwk8z.gt.1) ii8z = iiwk8z if (irwk8z.gt.1) ir8z = irwk8z C *******allocate workspace allocate (iwrk8z(ii8z),rwrk8z(ir8z)) plot = .false. C *******call pde solver call dtdp1x(xgrid, -1, neqn, nint, nbint, xout8z, uout, tout8z, ip &erdc, plot, lsqfit, fdiff, npts8z, t0, dt, nsteps, nout, nsave, cr &ankn, noupdt, itype, linear, rwrk8z, ir8z, iwrk8z, ii8z, restrt, g &ridid) deallocate (iwrk8z,rwrk8z) C *******call user-written postprocessor call postpr(tout8z,nsave,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 X (possibly as modified by UPRINT,..) # C 2 Y # C 3 U # C 4 V # C . . # C . . # C # C############################################################################## IVAR = 1 C T IS VARIABLE ics8z = 2 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, the plot will be scaled to just fit in the plot area. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## UMIN = -1.0 UMAX = 1.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 = 'X ' is8z = 0 ix8z = 0 call dtdpzp(ics8z,2*ivar-1,tout8z,nsave,xout8z,0,uout,neqn,title,u &min,umax,ix8z,is8z) C *******LINE PLOTS C############################################################################## C Enter a value for IVAR, to select the variable to be plotted or # C printed: # C IVAR = 1 means X (possibly as modified by UPRINT,..) # C 2 Y # C 3 U # C 4 V # C . . # C . . # C # C############################################################################## IVAR = 2 C T IS VARIABLE ics8z = 2 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, the plot will be scaled to just fit in the plot area. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## UMIN = -1.0 UMAX = 1.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 = 'Y ' is8z = 0 ix8z = 0 call dtdpzp(ics8z,2*ivar-1,tout8z,nsave,xout8z,0,uout,neqn,title,u &min,umax,ix8z,is8z) C *******LINE PLOTS C############################################################################## C Enter a value for IVAR, to select the variable to be plotted or # C printed: # C IVAR = 1 means X (possibly as modified by UPRINT,..) # C 2 Y # C 3 U # C 4 V # C . . # C . . # C # C############################################################################## IVAR = 3 C T IS VARIABLE ics8z = 2 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, the plot will be scaled to just fit in the plot area. +# 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 = 'U ' is8z = 0 ix8z = 0 call dtdpzp(ics8z,2*ivar-1,tout8z,nsave,xout8z,0,uout,neqn,title,u &min,umax,ix8z,is8z) C *******LINE PLOTS C############################################################################## C Enter a value for IVAR, to select the variable to be plotted or # C printed: # C IVAR = 1 means X (possibly as modified by UPRINT,..) # C 2 Y # C 3 U # C 4 V # C . . # C . . # C # C############################################################################## IVAR = 4 C T IS VARIABLE ics8z = 2 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, the plot will be scaled to just fit in the plot area. +# 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 = 'V ' is8z = 0 ix8z = 0 call dtdpzp(ics8z,2*ivar-1,tout8z,nsave,xout8z,0,uout,neqn,title,u &min,umax,ix8z,is8z) C *******LINE PLOTS C############################################################################## C Enter a value for IVAR, to select the variable to be plotted or # C printed: # C IVAR = 1 means X (possibly as modified by UPRINT,..) # C 2 Y # C 3 U # C 4 V # C . . # C . . # C # C############################################################################## IVAR = 5 C T IS VARIABLE ics8z = 2 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, the plot will be scaled to just fit in the plot area. +# 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 = 'W ' is8z = 0 ix8z = 0 call dtdpzp(ics8z,2*ivar-1,tout8z,nsave,xout8z,0,uout,neqn,title,u &min,umax,ix8z,is8z) 78755 continue call endgks stop end subroutine pdes8z(yd8z,i8z,j8z,kint8z,x8z,t,uu8z) 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 /dtdp4x/un8z(3,neqnmx) double precision uu8z(3,neqnmx) common/parm8z/ pi,RL ,G ,Rm X = uu8z(1, 1) Y = uu8z(1, 2) U = uu8z(1, 3) V = uu8z(1, 4) W = uu8z(1, 5) if (i8z.eq.0) then yd8z = 0.0 C############################################################################## C Enter FORTRAN expressions for the functions whose "integrals" are to # C be calculated and printed. They may be functions of X,Y,U,V # C ... and (if applicable) T. For 0D problems, the "integral" is just # C the value of the integrand. # C # C############################################################################## C 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 T,X,Y,U,V,... # C # C Recall that the PDEs have the form # C # C C11*d(X)/dT + C12*d(Y)/dT + C13*d(U)/dT + C14*d(V)/dT +...= F1 # C C21*d(X)/dT + C22*d(Y)/dT + C23*d(U)/dT + C24*d(V)/dT +...= F2 # C C31*d(X)/dT + C32*d(Y)/dT + C33*d(U)/dT + C34*d(V)/dT +...= F3 # C C41*d(X)/dT + C42*d(Y)/dT + C43*d(U)/dT + C44*d(V)/dT +...= F4 # C . . . . # C # C############################################################################## if (j8z.eq.0) then yd8z = 0.0 C-----------------------------------------> INPUT FROM GUI <------------------- C C(1,1) DEFINED if (i8z.eq. -101) yd8z = & 1 C C(1,2) DEFINED if (i8z.eq. -102) yd8z = & 0 C C(1,3) DEFINED if (i8z.eq. -103) yd8z = & 0 C C(1,4) DEFINED if (i8z.eq. -104) yd8z = & 0 C C(1,5) DEFINED if (i8z.eq. -105) yd8z = & 0 C-----------------------------------------> INPUT FROM GUI <------------------- C F1 DEFINED if (i8z.eq. 1) yd8z = & U C-----------------------------------------> INPUT FROM GUI <------------------- C C(2,1) DEFINED if (i8z.eq. -201) yd8z = & 0 C C(2,2) DEFINED if (i8z.eq. -202) yd8z = & 1 C C(2,3) DEFINED if (i8z.eq. -203) yd8z = & 0 C C(2,4) DEFINED if (i8z.eq. -204) yd8z = & 0 C C(2,5) DEFINED if (i8z.eq. -205) yd8z = & 0 C-----------------------------------------> INPUT FROM GUI <------------------- C F2 DEFINED if (i8z.eq. 2) yd8z = & V C-----------------------------------------> INPUT FROM GUI <------------------- C C(3,1) DEFINED if (i8z.eq. -301) yd8z = & 0 C C(3,2) DEFINED if (i8z.eq. -302) yd8z = & 0 C C(3,3) DEFINED if (i8z.eq. -303) yd8z = & Rm C C(3,4) DEFINED if (i8z.eq. -304) yd8z = & 0 C C(3,5) DEFINED if (i8z.eq. -305) yd8z = & 0 C-----------------------------------------> INPUT FROM GUI <------------------- C F3 DEFINED if (i8z.eq. 3) yd8z = & -W*X C-----------------------------------------> INPUT FROM GUI <------------------- C C(4,1) DEFINED if (i8z.eq. -401) yd8z = & 0 C C(4,2) DEFINED if (i8z.eq. -402) yd8z = & 0 C C(4,3) DEFINED if (i8z.eq. -403) yd8z = & 0 C C(4,4) DEFINED if (i8z.eq. -404) yd8z = & Rm C C(4,5) DEFINED if (i8z.eq. -405) yd8z = & 0 C-----------------------------------------> INPUT FROM GUI <------------------- C F4 DEFINED if (i8z.eq. 4) yd8z = & -W*Y - Rm*G C-----------------------------------------> INPUT FROM GUI <------------------- C C(5,1) DEFINED if (i8z.eq. -501) yd8z = & 0 C C(5,2) DEFINED if (i8z.eq. -502) yd8z = & 0 C C(5,3) DEFINED if (i8z.eq. -503) yd8z = & 0 C C(5,4) DEFINED if (i8z.eq. -504) yd8z = & 0 C C(5,5) DEFINED if (i8z.eq. -505) yd8z = & 0 C-----------------------------------------> INPUT FROM GUI <------------------- C F5 DEFINED if (i8z.eq. 5) yd8z = & X**2+Y**2-RL**2 else endif endif return end function u8z(i8z,x8z,t0) implicit double precision (a-h,o-z) common/parm8z/ pi,RL ,G ,Rm u8z = 0.0 C############################################################################## C Now the initial values must be defined using FORTRAN expressions. # C They may reference the initial time T0. # C############################################################################## C X0 DEFINED if (i8z.eq. 1) u8z = & RL C Y0 DEFINED if (i8z.eq. 2) u8z = & 0 C U0 DEFINED if (i8z.eq. 3) u8z = & 0 C V0 DEFINED if (i8z.eq. 4) u8z = & 0 C W0 DEFINED if (i8z.eq. 5) u8z = & 0 return end subroutine pmod8z(x8z,t,uu8z,uprint,uxp8z) implicit double precision (a-h,o-z) dimension uu8z(3,*),uprint(*),uxp8z(*) common/dtdp14/sint(20),bint(20),slim8z(20),blim8z(20) common/parm8z/ pi,RL ,G ,Rm X = uu8z(1, 1) Y = uu8z(1, 2) U = uu8z(1, 3) V = uu8z(1, 4) W = uu8z(1, 5) 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 X,Y,U,V,... 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 X +# C + UPRINT(2) Y +# C + UPRINT(3) U +# C + UPRINT(4) V +# C + . . +# C + . . +# C + Each function may be a function of X,Y,U,V,... and (if +# C + applicable) T, and of the "integral" estimates SINT(1),... +# C + +# C + The default for each variable is no change, for example, UPRINT(1) +# C + defaults to X. Enter FORTRAN expressions for each of the +# C + following functions (or default). +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C DEFINE UPRINT(*) HERE: return end C dummy routines subroutine gb8z(gd8z,ifac8z,i8z,j8z,x,t,uu8z) 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 function fb8z(i8z,iarc8z,ktri,s,x,y,t) implicit double precision (a-h,o-z) fb8z = 0 return end function axis8z(i8z,x,y,z,ical8z) implicit double precision (a-h,o-z) axis8z = 0 return end subroutine tran8z(itrans,x,y,z) implicit double precision (a-h,o-z) return end subroutine postpr(tout,nsave,uout,neqn) implicit double precision (a-h,o-z) dimension tout(0:nsave) dimension uout(2,neqn,0:nsave) common/parm8z/ pi,RL ,G ,Rm common /dtdp27/ itask,npes,icomm common /dtdp46/ eps8z,cgtl8z,npmx8z,itype,near8z data lun,lud/0,47/ if (itask.gt.0) return C UOUT(1,IEQ,L) = U_IEQ C (possibly as modified by UPRINT,..) 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! write (lun,*) 'fid = fopen(''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!78754 format (i8) C! do 78756 l=0,nsave0 C! write (lud,78757) tout(l) C! do 78755 ieq=1,neqn C! write (lud,78757) uout(1,ieq,l) C!78755 continue C!78756 continue C!78757 format (e16.8) C ******* WRITE pde2d.m C! call mtdp0d(itype,lun) return end