C ************************** C * PDE2D 9.2 MAIN PROGRAM * C ************************** C *** 2D 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) parameter (nxgrid=0,nygrid=0,npgrid=0,nqgrid=0) C############################################################################## C NV0 = number of vertices in initial triangulation # C############################################################################## PARAMETER (NV0 = 6) C############################################################################## C NT0 = number of triangles in initial triangulation # C############################################################################## PARAMETER (NT0 = 5) 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=101,NYP8Z=101,KDEG8Z=1,NBPT8Z=51) C############################################################################## C The solution is normally saved on a NX+1 by NY+1 rectangular grid of # C points # C (XA + I*(XB-XA)/NX , YA + J*(YB-YA)/NY) # C I=0,...,NX, J=0,...,NY. Enter values for NX and NY. Suggested values # C are NX=NY=25. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you want to save the solution at an arbitrary user-specified +# C + set of points, set NY=0 and NX+1=number of points. In this case you +# C + can request tabular output of the solution, but you cannot make any +# C + solution plots. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## PARAMETER (NX = 21) PARAMETER (NY = 21) PARAMETER (NSAVE = 1) common/parm8z/ pi,rho ,Rmu dimension vxy(2,nv0+1),iabc(3,nt0+1),iarc(nt0+1),xgrid(nxgrid+1),y &grid(nygrid+1),ixarc(2),iyarc(2),pgrid(npgrid+1),qgrid(nqgrid+1),i &parc(2),iqarc(2),xbd8z(nbpt8z,nt0+4),ybd8z(nbpt8z,nt0+4),xout8z(0: &nx,0:ny),yout8z(0:nx,0:ny),inrg8z(0:nx,0:ny),xcross(100),ycross(10 &0),tout8z(0:nsave),uout8z(0:nx,0:ny,3*neqn,0:nsave),uout(0:nx,0:ny &,3,neqn,0:nsave),xres8z(nxp8z),yres8z(nyp8z),ures8z(neqn,nxp8z,nyp &8z),xd0(ndelmx),yd0(ndelmx) equivalence (uout,uout8z) allocatable iwrk8z(:),rwrk8z(:) C dimension iwrk8z(iiwk8z),rwrk8z(irwk8z) character*40 title logical plot,symm,fdiff,evcmpx,crankn,noupdt,adapt,nodist,fillin,e &con8z,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/dtdp22/ nxa8z,nya8z,ifgr8z,kd8z,nbp8z common/dtdp23/ work8z(nxp8z*nyp8z+6) common/dtdp30/ econ8z,ncon8z common/dtdp46/ eps8z,cgtl8z,npmx8z,itype common/dtdp63/ amin8z(3*neqnmx),amax8z(3*neqnmx) common/dtdp65/ intri,iotri pi = 4.0*atan(1.d0) nxa8z = nxp8z nya8z = nyp8z kd8z = kdeg8z nbp8z = nbpt8z 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############################################################################## rho = & 1.1 Rmu = & 0.1 C *******INITIAL TRIANGULATION OPTION INTRI = 3 C *******SET IOTRI = 1 TO DUMP FINAL TRIANGULATION TO FILE pde2d.tri IOTRI = 0 C############################################################################## C For a general region, an initial triangulation is constructed # C which generally consists of only as many triangles as needed to define # C the region and to satisfy the following rules (triangles adjacent # C to a curved boundary may be considered to have one curved edge): # C # C 1. The end points of each arc are included as vertices in the # C triangulation. # C # C 2. No vertex of any triangle may touch another in a point which is # C not a vertex of the other triangle. # C # C 3. No triangle may have all three vertices on the boundary. # C # C Now enter the number of vertices in the initial triangulation (NV0) # C and the vertices (VXY(1,i),VXY(2,i)), i=1,...,NV0, when prompted. # C # C The vertices may be numbered in any order, but that order will define # C the vertex numbers referred to in the next list. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If the elements of VXY, IABC and IARC are defaulted, these initial +# C + triangulation arrays will be read from the file 'pde2d.tri'; in this +# C + case the values of NV0 and NT0 entered below must be the same as +# C + those in the file. +# C + +# C + A file 'pde2d.tri' is normally created when the final triangulation +# C + from another program (cases INTRI=1,2 or 3) is dumped. Reset IOTRI +# C + to 1 in the other program, to cause the final triangulation to be +# C + dumped into pde2d.tri. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## call dtdpwx(vxy,2*nv0,0) C VXY(1,1) = & -1 VXY(2,1) = & -1 C VXY(1,2) = & 1 VXY(2,2) = & -1 C VXY(1,3) = & 1 VXY(2,3) = & 1 C VXY(1,4) = & 0 VXY(2,4) = & 0.2 C VXY(1,5) = & -1 VXY(2,5) = & 1 C VXY(1,6) = & 0 VXY(2,6) = & -0.5 C############################################################################## C Now enter the number of triangles in the initial triangulation (NT0), # C and for each triangle k, give the vertex numbers of vertices a,b,c: # C IABC(1,k) , IABC(2,k) , IABC(3,k) # C where the third vertex (c) is not on the boundary of R (or on a curved # C interface arc), and the number, # C IARC(k) # C of the arc cut off by the base, ab, of the triangle. Put IARC(k)=0 # C if the base of triangle k does not intersect any boundary (or curved # C interface) arc. Recall that negative arc numbers correspond to # C 'fixed' boundary conditions, and positive arc numbers ( < 1000) # C correspond to 'free' boundary conditions. # C # C The order of this list defines the initial triangle numbers. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Any of the partial differential equation coefficients or other +# C + functions of X and Y which you are asked to supply later may be +# C + defined as functions of a variable 'KTRI', which holds the number +# C + of the INITIAL triangle in which (X,Y) lies. This is useful when +# C + the region is a composite of different materials, so that some +# C + material parameters have different values in different subregions. +# C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++# C + If any interfaces between subregions are curved, the interface arcs +# C + may be assigned unique arc numbers of 1000 and above, and treated +# C + like boundary arcs in the initial triangulation definition. (You +# C + will not be allowed to specify boundary conditions on them, however.)+# C + The triangulation refinement will follow the interface arcs, so that +# C + no final triangles will straddle an interface. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C IABC(1,1) = & 1 IABC(2,1) = & 2 IABC(3,1) = & 6 IARC(1) = & -1 C IABC(1,2) = & 2 IABC(2,2) = & 3 IABC(3,2) = & 6 IARC(2) = & 1 C IABC(1,3) = & 3 IABC(2,3) = & 4 IABC(3,3) = & 6 IARC(3) = & 1 C IABC(1,4) = & 4 IABC(2,4) = & 5 IABC(3,4) = & 6 IARC(4) = & 1 C IABC(1,5) = & 5 IABC(2,5) = & 1 IABC(3,5) = & 6 IARC(5) = & 1 call dtdpu(vxy,nv0,iabc,iarc,nt0) C############################################################################## C How many triangles (NTF) are desired for the final triangulation? # C############################################################################## NTF = 1000 C############################################################################## C If you don't want to read the FINE PRINT, enter ISOLVE = 4. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + The following linear system solvers are available: +# C + +# C + 1. Band method +# C + The band solver uses a reverse Cuthill-McKee ordering. +# C + 2. Frontal method +# C + This is an out-of-core version of the band solver. +# C + 3. Jacobi bi-conjugate gradient method +# C + This is a preconditioned bi-conjugate gradient, or +# C + Lanczos, iterative method. (This solver is MPI- +# C + enhanced, if MPI is available.) If you want to +# C + override the default convergence tolerance, set a +# C + new relative tolerance CGTL8Z in the main program. +# C + 4. Sparse direct method +# C + This is based on Harwell Library routines MA27/MA37, +# C + developed by AEA Industrial Technology at Harwell +# C + Laboratory, Oxfordshire, OX11 0RA, United Kingdom +# C + (used by permission). +# C + 5. Local solver +# C + Choose this option ONLY if alternative linear system +# C + solvers have been installed locally. See subroutines +# C + (D)TD3M, (D)TD3N in file (d)subs.f for instructions +# C + on how to add local solvers. +# C + 6. MPI-based parallel band solver +# C + This is a parallel solver which runs efficiently on +# C + multiple processor machines, under MPI. It is a +# C + band solver, with the matrix distributed over the +# C + available processors. Choose this option ONLY if the +# C + solver has been activated locally. See subroutine +# C + (D)TD3O in file (d)subs.f for instructions on how to +# C + activate this solver and the MPI-enhancements to the +# C + conjugate gradient solver. +# C + +# C + Enter ISOLVE = 1,2,3,4,5 or 6 to select a linear system solver. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## ISOLVE = 4 C############################################################################## C Enter the element degree (1,2,3 or 4) desired. A suggested value is # C IDEG = 3. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + A negative value for IDEG can be entered, and elements of degree +# C + ABS(IDEG) will be used, with a lower order numerical integration +# C + scheme. This results in a slight increase in speed, but negative +# C + values of IDEG are normally not recommended. However, IDEG = -1 +# C + is the only choice which will produce a lumped (diagonal) mass +# C + matrix. +# C + +# C + The spatial discretization error is O(h**2), O(h**3), O(h**4) or +# C + O(h**5) when elements of degree 1,2,3 or 4, respectively, are used, +# C + where h is the maximum triangle diameter, even if the region is +# C + curved. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## IDEG = 3 C *******STEADY-STATE PROBLEM itype = 1 t0 = 0.0 dt = 1.0 crankn = .false. noupdt = .false. C############################################################################## C Give an upper limit on the number of Newton's method iterations # C (NSTEPS) to be allowed for this nonlinear problem. NSTEPS defaults # C to 15. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + The iteration will stop if convergence occurs before the upper +# C + limit has been reached. In the FORTRAN program created by the +# C + preprocessor, NCON8Z will be .TRUE. if Newton's method has converged.+# C + +# C + For highly non-linear problems you may want to construct a one- +# C + parameter family of problems using the variable T, such that for +# C + T=1 the problems is easy (e.g. linear) and for T > N (N is less +# C + NSTEPS), the problem reduces to the original highly nonlinear +# C + problem. For example, the nonlinear term(s) may be multiplied by +# C + MIN(1.0,(T-1.0)/N). +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## NSTEPS = 15 NSTEPS = & 15 C############################################################################## C PDE2D solves the system of equations: # C # C d/dX* A1(X,Y,PHI,PHIx,PHIy,W,Wx,Wy) # C + d/dY* B1(X,Y,PHI,PHIx,PHIy,W,Wx,Wy) # C = F1(X,Y,PHI,PHIx,PHIy,W,Wx,Wy) # C # C d/dX* A2(X,Y,PHI,PHIx,PHIy,W,Wx,Wy) # C + d/dY* B2(X,Y,PHI,PHIx,PHIy,W,Wx,Wy) # C = F2(X,Y,PHI,PHIx,PHIy,W,Wx,Wy) # C # C with 'fixed' boundary conditions: # C # C PHI = FB1(X,Y) # C W = FB2(X,Y) # C # C or 'free' boundary conditions: # C # C A1*nx + B1*ny = GB1(X,Y,PHI,PHIx,PHIy,W,Wx,Wy) # C A2*nx + B2*ny = GB2(X,Y,PHI,PHIx,PHIy,W,Wx,Wy) # C # C where PHI(X,Y) and W(X,Y) are the unknowns and F1,A1,B1,F2,A2,B2, # C FB1,FB2,GB1,GB2 are user-supplied functions. # C # C note: # C (nx,ny) = unit outward normal to the boundary # C PHIx = d(PHI)/dX Wx = d(W)/dX # C PHIy = d(PHI)/dY Wy = d(W)/dY # C # C Is this problem symmetric? If you don't want to read the FINE PRINT, # C it is safe to enter 'no'. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + This problem is called symmetric if each of the matrices +# C + +# C + F1.PHI F1.PHIx F1.PHIy F1.W F1.Wx F1.Wy +# C + A1.PHI A1.PHIx A1.PHIy A1.W A1.Wx A1.Wy +# C + B1.PHI B1.PHIx B1.PHIy B1.W B1.Wx B1.Wy +# C + F2.PHI F2.PHIx F2.PHIy F2.W F2.Wx F2.Wy +# C + A2.PHI A2.PHIx A2.PHIy A2.W A2.Wx A2.Wy +# C + B2.PHI B2.PHIx B2.PHIy B2.W B2.Wx B2.Wy +# C + +# C + and +# C + GB1.PHI GB1.W +# C + GB2.PHI GB2.W +# C + +# C + is always symmetric, where F1.PHI means d(F1)/d(PHI), and similarly +# C + for the other terms. In addition, GB1,GB2 must not depend on +# C + PHIx,PHIy,Wx,Wy. +# C + +# C + The memory and execution time are halved if the problem is known to +# C + be symmetric. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## SYMM = .FALSE. C############################################################################## C If you don't want to read the FINE PRINT, enter 'yes' (strongly # C recommended). # 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 = 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 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. C GRIDID = .FALSE. IF FINITE ELEMENT GRID CHANGES BETWEEN DUMP, RESTART GRIDID = .TRUE. C############################################################################## C The solution is saved on an NX+1 by NY+1 rectangular grid covering the # C rectangle (XA,XB) x (YA,YB). Enter values for XA,XB,YA,YB. These # C variables are usually defaulted. # C # C The default is a rectangle which just covers the entire region. # C############################################################################## C defaults for xa,xb,ya,yb call dtdpv (xgrid,nxgrid,ygrid,nygrid,vxy,nv0,iarc,nt0,xa,xb,ya, &yb) C DEFINE XA,XB,YA,YB HERE: call dtdpx2(nx,ny,xa,xb,ya,yb,hx8z,hy8z,xout8z,yout8z,npts8z) C SOLUTION SAVED EVERY NOUT ITERATIONS NOUT = NSTEPS C *******allocate workspace call dtdpzz(ntf,ideg,isolve,symm,neqn,ii8z,ir8z) if (iiwk8z.gt.1) ii8z = iiwk8z if (irwk8z.gt.1) ir8z = irwk8z allocate (iwrk8z(ii8z),rwrk8z(ir8z)) C *******DRAW TRIANGULATION PLOTS (OVER C *******RECTANGLE (XA,XB) x (YA,YB))? PLOT = .TRUE. C *******call pde solver call dtdp2x(xgrid, nxgrid, ygrid, nygrid, ixarc, iyarc, vxy, nv0, &iabc, nt0, iarc, restrt, gridid, neqn, ntf, ideg, isolve, nsteps, &nout, t0, dt, plot, symm, fdiff, itype, nint, nbint, ndel, xd0, yd &0, crankn, noupdt, xbd8z, ybd8z, nbd8z, xout8z, yout8z, uout, inrg &8z, npts8z, ny, tout8z, nsave, iwrk8z, ii8z, rwrk8z, ir8z) deallocate (iwrk8z,rwrk8z) C *******UNCOMMENT TO POSTPROCESS ONLY IF CONVERGED C IF (.NOT.NCON8Z) GO TO 78755 C *******read from restart file to array ures8z C call dtdpr2(1,xres8z,nxp8z,yres8z,nyp8z,ures8z,neqn) C *******write array ures8z back to restart file C call dtdpr2(2,xres8z,nxp8z,yres8z,nyp8z,ures8z,neqn) C *******call user-written postprocessor call postpr(tout8z,nsave,xout8z,yout8z,inrg8z,nx,ny,uout,neqn) C *******VECTOR FIELD PLOT C############################################################################## C Enter values for IVARX, IVARY to select the X and Y components of # C the vector to be plotted. # C IVARX or IVARY = 1 means PHI (possibly as modified by UPRINT,..) # C 2 A1 # C 3 B1 # C 4 W # C 5 A2 # C 6 B2 # C############################################################################## IVARX = 2 IVARY = 3 ISET1 = 1 ISET2 = NSAVE ISINC = 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 + Do you want to scale the axes on the plot so that the region is +# C + undistorted? Otherwise the axes will be scaled so that the figure +# C + approximately fills the plot space. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## NODIST = .TRUE. C a1mag = max(abs(amin8z(ivarx)),abs(amax8z(ivarx))) a2mag = max(abs(amin8z(ivary)),abs(amax8z(ivary))) C############################################################################## C For the purpose of scaling the arrows, the ranges of the two components # C of the vector are assumed to be (-VR1MAG,VR1MAG) and (-VR2MAG,VR2MAG). # C Enter values for VR1MAG and VR2MAG. VR1MAG and VR2MAG are often # C defaulted. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + By default, VR1MAG and VR2MAG are the maxima of the absolute values +# C + of the first and second components. For a common scaling, you may +# C + want to set VR1MAG=A1MAG, VR2MAG=A2MAG. A1MAG, A2MAG are the +# C + maxima of the absolute values over all output points and over all +# C + saved time steps or iterations. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## VR1MAG = 0.0 VR2MAG = 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 = 'Fluid Velocity (PHIy,-PHIx) ' call dtdprx(tout8z,nsave,iset1,iset2,isinc) do 78756 is8z=iset1,iset2,isinc call dtdpla(uout8z(0,0,ivarx,is8z),uout8z(0,0,ivary,is8z),nx,ny,xa &,ya,hx8z,hy8z,inrg8z,xbd8z,ybd8z,nbd8z,title,vr1mag,vr2mag,nodist, &tout8z(is8z)) 78756 continue C *******CONTOUR PLOT C############################################################################## C Enter a value for IVAR, to select the variable to be plotted or # C printed: # C IVAR = 1 means PHI (possibly as modified by UPRINT,..) # C 2 A1 # C 3 B1 # C 4 W # C 5 A2 # C 6 B2 # C############################################################################## IVAR = 1 ISET1 = 1 ISET2 = NSAVE ISINC = 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 + Do you want to scale the axes on the plot so that the region is +# C + undistorted? Otherwise the axes will be scaled so that the figure +# C + approximately fills the plot space. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## NODIST = .TRUE. C alow = amin8z(ivar) ahigh = amax8z(ivar) C############################################################################## C Enter lower (UMIN) and upper (UMAX) bounds for the contour values. UMIN # C and UMAX are often defaulted. # C # C Labeled contours will be drawn corresponding to the values # C # C UMIN + S*(UMAX-UMIN), for S=0.05,0.15,...0.95. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + By default, UMIN and UMAX are set to the minimum and maximum values +# C + of the variable to be plotted. For a common scaling, you may want +# C + to set UMIN=ALOW, UMAX=AHIGH. ALOW and AHIGH are the minimum and +# C + maximum values over all output points and over all saved time steps +# C + or iterations. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## UMIN = 0.0 UMAX = 0.0 C############################################################################## C Do you want two additional unlabeled contours to be drawn between each # C pair of labeled contours? # C############################################################################## FILLIN = .FALSE. C############################################################################## C Enter a title, WITHOUT quotation marks. A maximum of 40 characters # C are allowed. The default is no title. # C############################################################################## TITLE = ' ' TITLE = 'Stream function (PHI) ' call dtdprx(tout8z,nsave,iset1,iset2,isinc) do 78757 is8z=iset1,iset2,isinc call dtdplg(uout8z(0,0,ivar,is8z),nx,ny,xa,ya,hx8z,hy8z,inrg8z,xbd &8z,ybd8z,nbd8z,title,umin,umax,nodist,fillin,tout8z(is8z)) 78757 continue C *******CONTOUR PLOT C############################################################################## C Enter a value for IVAR, to select the variable to be plotted or # C printed: # C IVAR = 1 means PHI (possibly as modified by UPRINT,..) # C 2 A1 # C 3 B1 # C 4 W # C 5 A2 # C 6 B2 # C############################################################################## IVAR = 4 ISET1 = 1 ISET2 = NSAVE ISINC = 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 + Do you want to scale the axes on the plot so that the region is +# C + undistorted? Otherwise the axes will be scaled so that the figure +# C + approximately fills the plot space. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## NODIST = .TRUE. C alow = amin8z(ivar) ahigh = amax8z(ivar) C############################################################################## C Enter lower (UMIN) and upper (UMAX) bounds for the contour values. UMIN # C and UMAX are often defaulted. # C # C Labeled contours will be drawn corresponding to the values # C # C UMIN + S*(UMAX-UMIN), for S=0.05,0.15,...0.95. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + By default, UMIN and UMAX are set to the minimum and maximum values +# C + of the variable to be plotted. For a common scaling, you may want +# C + to set UMIN=ALOW, UMAX=AHIGH. ALOW and AHIGH are the minimum and +# C + maximum values over all output points and over all saved time steps +# C + or iterations. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## UMIN = 0.0 UMAX = 0.0 C############################################################################## C Do you want two additional unlabeled contours to be drawn between each # C pair of labeled contours? # C############################################################################## FILLIN = .FALSE. C############################################################################## C Enter a title, WITHOUT quotation marks. A maximum of 40 characters # C are allowed. The default is no title. # C############################################################################## TITLE = ' ' TITLE = 'Vorticity (W) ' c call dtdprx(tout8z,nsave,iset1,iset2,isinc) c do 78758 is8z=iset1,iset2,isinc c call dtdplg(uout8z(0,0,ivar,is8z),nx,ny,xa,ya,hx8z,hy8z,inrg8z,xbd c &8z,ybd8z,nbd8z,title,umin,umax,nodist,fillin,tout8z(is8z)) c78758 continue 78755 continue call endgks stop end subroutine tran8z(p,q,x,y) implicit double precision (a-h,o-z) x = p y = q return end subroutine xy8z(i8z,iarc8z,s,x,y,s0,sf) implicit double precision (a-h,o-z) dimension pxy(2,1000) common/parm8z/ pi,rho ,Rmu x = 0.0 y = 0.0 return end subroutine dis8z(x,y,ktri,triden,shape) implicit double precision (a-h,o-z) common/parm8z/ pi,rho ,Rmu adapt = dtdplx(2) C############################################################################## C Now enter a FORTRAN expression for TRIDEN(X,Y), which controls the # C grading of the triangulation. TRIDEN should be largest where the # C triangulation is to be most dense. The default is TRIDEN(X,Y)=1.0 # C (a uniform triangulation). # C # C If you want the triangulation to be graded "adaptively", set # C TRIDEN = ADAPT and and make sure there is no "pde2d.adp" file in the # C working directory the first time you run the program, then run the # C program 2 or more times, possibly increasing NTF each time. On the # C first run, a uniform triangulation will be generated, but on each # C subsequent run, information output by the previous run to "pde2d.adp" # C will be used to guide the grading of the new triangulation. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If TRIDEN = ADAPT, after each run a file "pde2d.adp" is written +# C + which tabulates the values of the magnitude of the gradient of the +# C + solution (at the last time step or iteration) at an output NXP8Z by +# C + NYP8Z grid of points (NXP8Z and NYP8Z are set to 101 in a PARAMETER +# C + statement in the main program, so they can be changed if desired). +# C + If NEQN > 1, a normalized average of the gradients of the NEQN +# C + solution components is used. +# C + +# C + You can do all the "runs" in one program, by setting NPROB > 1. +# C + Each pass through the DO loop, PDE2D will read the gradient values +# C + output the previous pass. +# C + +# C + TRIDEN may also be a function of the initial triangle number KTRI. +# C + +# C + Increase the variable EXAG from its default value of 1.5 if you want +# C + to exaggerate the grading of an adaptive triangulation (make it less +# C + uniform). EXAG should normally not be larger than about 2.0. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## TRIDEN = 1.0 EXAG = 1.5 C############################################################################## C If you don't want to read the FINE PRINT, default SHAPE. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Enter a FORTRAN expression for SHAPE(X,Y), which controls the +# C + approximate shape of the triangles. The triangulation refinement +# C + will proceed with the goal of generating triangles with an average +# C + height to width ratio of approximately SHAPE(X,Y) near the point +# C + (X,Y). SHAPE must be positive. The default is SHAPE(X,Y)=1.0. +# C + +# C + SHAPE may also be a function of the initial triangle number KTRI. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## SHAPE = 1.0 if (triden.eq.adapt) triden = dtdpgr(x,y)**exag return end subroutine pdes8z(yd8z,i8z,j8z,kint8z,idel8z,jdel8z,x,y,ktri,t) implicit double precision (a-h,o-z) parameter (neqnmx= 99) parameter (ndelmx= 20) C un8z(1,I),un8z(2,I),un8z(3,I) hold the (rarely used) values C of UI,UIx,UIy from the previous iteration or time step common/dtdp4/un8z(3,neqnmx),uu8z(3,neqnmx) common/dtdp17/normx,normy,iarc double precision normx,normy,delamp(ndelmx,neqnmx) common/parm8z/ pi,rho ,Rmu PHI = uu8z(1, 1) PHIx = uu8z(2, 1) PHIy = uu8z(3, 1) PHInorm = PHIx*normx + PHIy*normy W = uu8z(1, 2) Wx = uu8z(2, 2) Wy = uu8z(3, 2) Wnorm = Wx*normx + Wy*normy 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,Y,PHI,PHIx,PHIy,W,Wx,Wy and (if applicable) T # C # C The integrals may also contain references to the initial triangle # C number KTRI. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you only want to integrate a function over part of the region, +# C + define that function to be zero in the rest of the region. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C INTEGRAL DEFINED if (kint8z.eq. 1) yd8z = & W 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,Y,PHI,PHIx,PHIy,W,Wx,Wy and (if applicable) T # C # C The components (NORMx,NORMy) of the unit outward normal vector, and the # C initial triangle number KTRI, and the boundary arc number IARC may also # C be referenced. You can also reference the normal derivatives PHInorm, # C Wnorm. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + CURVED interior interface arcs are considered part of the boundary, +# C + for the boundary integral computations, ONLY IF they have arc numbers+# C + in the range 8000-8999. In this case, since an interface arc is +# C + considered to be a boundary for both of the subregions it separates, +# C + the boundary integral will be computed twice on each curved +# C + interface arc, once with (NORMx,NORMy) defined in each direction. +# C + +# C + If you only want to integrate a function over part of the boundary, +# C + define that function to be zero on the rest of the boundary. You +# C + can examine the point (X,Y) to determine if it is on the desired +# C + boundary segment, or the boundary arc number IARC, or the initial +# C + triangle number KTRI (if INTRI=3), +# 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,Y,PHI,PHIx,PHIy,W,Wx,Wy # C # C They may also be functions of the initial triangle number KTRI # C and, in some cases, of the parameter T. # C # C Recall that the PDEs have the form # C # C d/dX*A1 + d/dY*B1 = F1 # C d/dX*A2 + d/dY*B2 = F2 # C # C############################################################################## F1y = 1 F2x = -1 beta = min(1.d0,(T-1)/5.0) if (j8z.eq.0) then yd8z = 0.0 C F1 DEFINED if (i8z.eq. 1) yd8z = & W C A1 DEFINED if (i8z.eq. 2) yd8z = & PHIx C B1 DEFINED if (i8z.eq. 3) yd8z = & PHIy C F2 DEFINED if (i8z.eq. 4) yd8z = & beta*rho*(PHIy*Wx-PHIx*Wy) + F2x - F1y C A2 DEFINED if (i8z.eq. 5) yd8z = & Rmu*Wx C B2 DEFINED if (i8z.eq. 6) yd8z = & Rmu*Wy else endif endif return end function u8z(i8z,x,y,ktri,t0) implicit double precision (a-h,o-z) common/parm8z/ pi,rho ,Rmu u8z = 0.0 C############################################################################## C Now the initial values for Newton's method must be defined using # C FORTRAN expressions. They may be functions of X and Y, and of the # C initial triangle number KTRI. # C # C It is important to provide initial values which are at least of the # C correct order of magnitude. # C############################################################################## C PHI0 DEFINED if (i8z.eq. 1) u8z = & 0 C W0 DEFINED if (i8z.eq. 2) u8z = & 0 return end function fb8z(i8z,iarc8z,ktri,s,x,y,t) implicit double precision (a-h,o-z) common/parm8z/ pi,rho ,Rmu fb8z = 0.0 C############################################################################## C Enter the arc number (IARC) of a boundary arc with negative arc number. # C############################################################################## IARC = -1 if (iarc8z.eq.iarc) then C############################################################################## C Enter FORTRAN expressions to define FB1,FB2 on this arc. They may # C be functions of X,Y and (if applicable) T. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + These functions may also reference the initial triangle number KTRI, +# C + and the arc parameter S (S = P or Q when INTRI=2). +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C # C Recall that fixed boundary conditions have the form # C # C PHI = FB1 # C W = FB2 # C # C############################################################################## C FB1 DEFINED if (i8z.eq. 1) fb8z = & 0 C FB2 DEFINED if (i8z.eq. 2) fb8z = & 0 endif return end subroutine gb8z(gd8z,i8z,j8z,iarc8z,ktri,s,x,y,t) implicit double precision (a-h,o-z) parameter (neqnmx= 99) C un8z(1,I),un8z(2,I),un8z(3,I) hold the (rarely used) values C of UI,UIx,UIy from the previous iteration or time step. C (normx,normy) is the (rarely used) unit outward normal vector common/dtdp4/un8z(3,neqnmx),uu8z(neqnmx,3) common/dtdp17/normx,normy,ibarc8z common/dtdp49/bign8z double precision normx,normy common/parm8z/ pi,rho ,Rmu zero(f8z) = bign8z*f8z PHI = uu8z( 1,1) PHIx = uu8z( 1,2) PHIy = uu8z( 1,3) W = uu8z( 2,1) Wx = uu8z( 2,2) Wy = uu8z( 2,3) if (j8z.eq.0) gd8z = 0.0 C############################################################################## C Enter the arc number (IARC) of a boundary arc with positive arc number. # C############################################################################## IARC = 1 if (iarc8z.eq.iarc) then C############################################################################## C Enter FORTRAN expressions to define the following free boundary # C condition functions on this arc. They may be functions of # C # C X,Y,PHI,PHIx,PHIy,W,Wx,Wy and (if applicable) T # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + These functions may also reference the components (NORMx,NORMy) of +# C + the unit outward normal vector, the initial triangle number KTRI, +# C + and the arc parameter S (S = P or Q when INTRI=2). +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C # C Recall that free boundary conditions have the form # C # C A1*nx+B1*ny = GB1 # C A2*nx+B2*ny = GB2 # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Note that 'fixed' boundary conditions: +# C + Ui = FBi(X,Y,[T]) +# C + can be expressed as 'free' boundary conditions in the form: +# C + GBi = zero(Ui-FBi(X,Y,[T])) +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## if (j8z.eq.0) then C GB1 DEFINED if (i8z.eq. 1) gd8z = & 0 C GB2 DEFINED if (i8z.eq. 2) gd8z = & zero(PHI) else endif endif return end subroutine pmod8z(x,y,ktri,t,a,b) implicit double precision (a-h,o-z) parameter (neqnmx= 99) common/dtdp4/un8z(3,neqnmx),uu8z(3,neqnmx) common/dtdp6/upr8z(3,neqnmx),uab8z(3,neqnmx) common/dtdp9/uprint(neqnmx),aprint(neqnmx),bprint(neqnmx) common/dtdp14/sint(20),bint(20),slim8z(20),blim8z(20) common/parm8z/ pi,rho ,Rmu PHI = uu8z(1, 1) PHIx = uu8z(2, 1) PHIy = uu8z(3, 1) a1 = upr8z(2, 1) b1 = upr8z(3, 1) W = uu8z(1, 2) Wx = uu8z(2, 2) Wy = uu8z(3, 2) a2 = upr8z(2, 2) b2 = upr8z(3, 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 PHI,A1,B1,W,A2,B2 at the +# C + output points. If different variables are to be saved (for later +# C + printing or plotting) the following functions can be used to +# C + re-define the output variables: +# C + define UPRINT(1) to replace PHI +# C + APRINT(1) A1 +# C + BPRINT(1) B1 +# C + UPRINT(2) W +# C + APRINT(2) A2 +# C + BPRINT(2) B2 +# C + Each function may be a function of +# C + +# C + X,Y,PHI,PHIx,PHIy,A1,B1,W,Wx,Wy,A2,B2 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 PHI. Enter FORTRAN expressions for each of the +# C + following functions (or default). +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C DEFINE UPRINT(*),APRINT(*),BPRINT(*) HERE: APRINT(1) = & PHIy BPRINT(1) = & -PHIx return end C dummy routines function axis8z(i8z,x,y,z,ical8z) implicit double precision (a-h,o-z) axis8z = 0 return end subroutine postpr(tout,nsave,xout,yout,inrg,nx,ny,uout,neqn) implicit double precision (a-h,o-z) dimension xout(0:nx,0:ny),yout(0:nx,0:ny),tout(0:nsave) dimension inrg(0:nx,0:ny),uout(0:nx,0:ny,3,neqn,0:nsave) common/parm8z/ pi,rho ,Rmu common /dtdp27/ itask,npes,icomm common /dtdp46/ eps8z,cgtl8z,npmx8z,itype data lun,lud/0,47/ if (itask.gt.0) return C UOUT(I,J,IDER,IEQ,L) = U-sub-IEQ, if IDER=1 C A-sub-IEQ, if IDER=2 C B-sub-IEQ, if IDER=3 C (possibly as modified by UPRINT,..) C at the point (XOUT(I,J) , YOUT(I,J)) C at time/iteration TOUT(L). C INRG(I,J) = 1 if this point is in R C = 0 otherwise 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! write (lud,78754) ny C!78754 format (i8) C! do 78756 i=0,nx C! do 78755 j=0,ny C! write (lud,78762) xout(i,j),yout(i,j) C!78755 continue C!78756 continue C! do 78761 l=0,nsave0 C! write (lud,78762) tout(l) C! do 78760 ieq=1,neqn C! do 78759 ider=1,3 C! do 78758 i=0,nx C! do 78757 j=0,ny C! if (inrg(i,j).eq.1) then C! write (lud,78762) uout(i,j,ider,ieq,l) C! else C! write (lud,78763) C! endif C!78757 continue C!78758 continue C!78759 continue C!78760 continue C!78761 continue C!78762 format (e16.8) C!78763 format ('NaN') 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! write (lun,*) 'NY = 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(NY+1,NX+1);' C! write (lun,*) 'Y = zeros(NY+1,NX+1);' C! write (lun,*) 'U = zeros(NY+1,NX+1,NSAVE+1,3,NEQN);' C! write (lun,*) 'for i=0:NX' C! write (lun,*) 'for j=0:NY' C! write (lun,*) ' X(j+1,i+1) = fscanf(fid,''%g'',1);' C! write (lun,*) ' Y(j+1,i+1) = fscanf(fid,''%g'',1);' C! write (lun,*) 'end' 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:3' C! write (lun,*) 'for i=0:NX' C! write (lun,*) 'for j=0:NY' C! write (lun,*) C! & ' U(j+1,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,*) 'end' C! write (lun,*) 'xmin = min(min(X(:,:)));' C! write (lun,*) 'xmax = max(max(X(:,:)));' C! write (lun,*) 'ymin = min(min(Y(:,:)));' C! write (lun,*) 'ymax = max(max(Y(:,:)));' C! write (lun,*) 'hx = 0.1*(xmax-xmin);' C! write (lun,*) 'hy = 0.1*(ymax-ymin);' C! write (lun,*) '% Surface plots of each variable' C! write (lun,*) 'for IEQ=1:NEQN' C! write (lun,*) 'IDER = 1;' C! write (lun,*) C! & 'umin = min(min(min(U(:,:,L0+1:NSAVE+1,IDER,IEQ))));' C! write (lun,*) C! & 'umax = max(max(max(U(:,:,L0+1:NSAVE+1,IDER,IEQ))));' C! write (lun,*) 'for L=L0:NSAVE' C! write (lun,*) ' figure' C! write (lun,*) ' surf(X,Y,U(:,:,L+1,IDER,IEQ))' C! write (lun,*) ' axis([xmin xmax ymin ymax umin umax])' C! write (lun,*) ' xlabel(''X'')' C! write (lun,*) ' ylabel(''Y'')' C! write (lun,*) ' zlabel([''U'',num2str(IEQ)])' C! write (lun,*) ' title(['' T = '',num2str(T(L+1))])' C! write (lun,*) ' view(-37.5,30.0)' C! write (lun,*) 'end' C! write (lun,*) 'end' C! write (lun,*) '% ******* Choose variables for vector plots' C! write (lun,*) '% (see comments in POSTPR)' C! write (lun,*) 'IEQ1 = 1;' C! write (lun,*) 'IDER1 = 2;' C! write (lun,*) 'IEQ2 = 1;' C! write (lun,*) 'IDER2 = 3;' C! write (lun,*) 'for L=L0:NSAVE' C! write (lun,*) ' figure' C! write (lun,*) ' quiver(X,Y,U(:,:,L+1,IDER1,IEQ1), ...' C! write (lun,*) ' U(:,:,L+1,IDER2,IEQ2))' C! write (lun,*) ' axis([xmin-hx xmax+hx ymin-hy ymax+hy])' C! write (lun,*) ' xlabel(''X'')' C! write (lun,*) ' ylabel(''Y'')' C! write (lun,*) ' title(['' T = '',num2str(T(L+1))])' C! write (lun,*) 'end' return end