C ************************** C * PDE2D 9.4 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 (nv0=0,nt0=0,npgrid=0,nqgrid=0) C############################################################################## C NXGRID = number of X-grid lines # C############################################################################## PARAMETER (NXGRID = 31) C############################################################################## C NYGRID = number of Y-grid lines # C############################################################################## PARAMETER (NYGRID = 21) C############################################################################## C How many differential equations (NEQN) are there in your problem? # C############################################################################## PARAMETER (NEQN = 3) C DIMENSIONS OF WORK ARRAYS C SET TO 1 FOR AUTOMATIC ALLOCATION PARAMETER (IRWK8Z= 50000000) PARAMETER (IIWK8Z= 1) PARAMETER (NXP8Z=1001,NYP8Z=1001,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 + +# C + If you set NEAR8Z=1 in the main program, the values saved at each +# C + output point will actually be the solution as evaluated at a nearby +# C + integration point. For most problems this obviously will produce +# C + less accurate output or plots, but for certain (rare) problems, a +# C + solution component may be much less noisy when plotted only at +# C + integration points. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## PARAMETER (NX = 80) PARAMETER (NY = 60) 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 = 1) common /parm8z/ pi,r1i,r1r,r2i,r2r,x01,y01,x02,y02 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),xd0(ndelmx),yd0(ndelmx) C dimension xres8z(nxp8z),yres8z(nyp8z),ures8z(neqn,nxp8z,nyp8z) 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,near8z common/dtdp63/ amin8z(3*neqnmx),amax8z(3*neqnmx) common/dtdp65/ intri,iotri common/dtdp76/ mdim8z,nx18z,ny18z,xa,xb,ya,yb,uout(0:nx,0:ny,3,neq &n,0:nsave) pi = 4.0*atan(1.d0) nxa8z = nxp8z nya8z = nyp8z nx18z = nx+1 ny18z = ny+1 mdim8z = 3 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 = 40 DDT = 0.125 do 78755 iprob=1,nprob C############################################################################## C PDE2D solves the time-dependent system (note: U,A,B,F,FB,GB,U0 may be # C vectors, C,RHO may be matrices): # C # C C(X,Y,T,U,Ux,Uy)*d(U)/dT = d/dX* A(X,Y,T,U,Ux,Uy) # C + d/dY* B(X,Y,T,U,Ux,Uy) # C - F(X,Y,T,U,Ux,Uy) # C # C or the steady-state system: # C # C d/dX* A(X,Y,U,Ux,Uy) # C + d/dY* B(X,Y,U,Ux,Uy) # C = F(X,Y,U,Ux,Uy) # C # C or the linear and homogeneous eigenvalue system: # C # C d/dX* A(X,Y,U,Ux,Uy) # C + d/dY* B(X,Y,U,Ux,Uy) # C = F(X,Y,U,Ux,Uy) + lambda*RHO(X,Y)*U # C # C in an arbitrary two-dimensional region, R, with 'fixed' boundary # C conditions on part of the boundary: # C # C U = FB(X,Y,[T]) # C # C and 'free' boundary conditions on the other part: # C # C A*nx + B*ny = GB(X,Y,[T],U,Ux,Uy) # C # C For time-dependent problems there are also initial conditions: # C # C U = U0(X,Y) at T=T0 # C # C Here Ux,Uy represent the (vector) functions dU/dX,dU/dY, and (nx,ny) # C represents the unit outward normal to the boundary. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + A system of NEQN complex partial differential equations must be +# C + written as a system of 2*NEQN real equations, by separating the +# C + equations into their real and imaginary parts. However, note that +# C + the complex arithmetic abilities of FORTRAN can be used to simplify +# C + this separation. For example, the complex PDE: +# C + I*(Uxx+Uyy) = 1/(1+U**10), where U = UR + UI*I +# C + would be difficult to split up analytically, but using FORTRAN +# C + expressions it is easy: +# C + A1 = -UIx, B1 = -UIy, F1 = REAL(1.0/(1.0+CMPLX(UR,UI)**10)) +# C + A2 = URx, B2 = URy, F2 = AIMAG(1.0/(1.0+CMPLX(UR,UI)**10)) +# C + +# C + If your PDEs involve the solution at points other than (X,Y), the +# C + function +# C + (D)OLDSOL2(IDER,IEQ,XX,YY,KDEG) +# C + will interpolate (using interpolation of degree KDEG=1,2 or 3) to +# C + (XX,YY) the function saved in UOUT(*,*,IDER,IEQ,ISET) on the last +# C + time step or iteration (ISET) for which it has been saved. Thus, +# C + for example, if IDER=1, this will return the latest value of +# C + component IEQ of the solution at (XX,YY), assuming this has not been +# C + modified using UPRINT... If your equations involve integrals of the +# C + solution, for example, you can use (D)OLDSOL2 to approximate these +# C + using the solution from the last time step or iteration. +# C + +# C + CAUTION: For a steady-state or eigenvalue problem, you must reset +# C + NOUT=1 if you want to save the solution each iteration. +# 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############################################################################## r1i = 0.5 r1r = r1i r2i = -0.5 r2r = r2i x01 = 15 y01 = -15 x02 = 15 y02 = 15 C *******INITIAL TRIANGULATION OPTION INTRI = 1 C *******SET IOTRI = 1 TO DUMP FINAL TRIANGULATION TO FILE pde2d.tri IOTRI = 0 C############################################################################## C For a rectangular region, the initial triangulation is defined by a # C set of grid lines corresponding to # C X = XGRID(1),XGRID(2),...,XGRID(NXGRID) # C Y = YGRID(1),YGRID(2),...,YGRID(NYGRID) # C Each grid rectangle is divided into four equal area triangles. # C # 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 represent the values of X on the left and right sides of the # C rectangle R, and cannot be defaulted. Then you will be prompted # C similarly for the number and values of the Y-grid points. # C############################################################################## call dtdpwx(xgrid,nxgrid,0) call dtdpwx(ygrid,nygrid,0) C XGRID DEFINED XGRID(1) = & 0 XGRID(NXGRID) = & 70 C YGRID DEFINED YGRID(1) = & -30 YGRID(NYGRID) = & 30 call dtdpwx(xgrid,nxgrid,1) call dtdpwx(ygrid,nygrid,1) C############################################################################## C Enter the arc numbers IXARC(1),IXARC(2) of the left and right sides, # C respectively, and the arc numbers IYARC(1),IYARC(2) of the bottom and # C top sides of the rectangle R. Recall that negative arc numbers # C correspond to 'fixed' boundary conditions, and positive arc numbers # C ( < 1000) correspond to 'free' boundary conditions. # C############################################################################## IXARC(1) = & -2 IXARC(2) = & 1 IYARC(1) = & 2 IYARC(2) = & 2 C############################################################################## C How many triangles (NTF) are desired for the final triangulation? # C############################################################################## NTF = 8*(NXGRID-1)*(NYGRID-1) 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. +# C + +# C + The spatial discretization error is O(h**2), O(h**3), O(h**4) or +# C + O(h**5) when IDEG = 1,2,3 or 4, respectively, is used, where h is +# C + the maximum triangle diameter, even if the region is curved. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## IDEG = 3 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 T0 = (iprob-1)*DDT TF = iprob*DDT 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############################################################################## ADAPT = .TRUE. TOLER(1) = 1.D10 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############################################################################## NSTEPS = 10 dt = (tf-t0)/max(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 the Crank-Nicolson scheme to be used to discretize time? If you +# C + answer 'no', a backward Euler scheme will be used. +# C + +# C + If a user-specified constant time step is chosen, the second order +# C + Crank Nicolson method is recommended only for problems with very +# C + well-behaved solutions, and the first order backward Euler scheme +# C + should be used for more difficult problems. In particular, do not +# C + use the Crank Nicolson method if the left hand side of any PDE is +# C + zero, for example, if a mixed elliptic/parabolic problem is solved. +# C + +# C + If adaptive time step control is chosen, however, an extrapolation +# C + is done between the 1-step and 2-step answers which makes the Euler +# C + method second order, and the Crank-Nicolson method strongly stable. +# C + Thus in this case, both methods have second order accuracy, and both +# C + are strongly stable. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## CRANKN = .TRUE. C############################################################################## C PDE2D solves the system of equations: # C # C C11(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy)*d(U)/dT # C + C12(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy)*d(V)/dT # C + C13(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy)*d(W)/dT = # C d/dX* A1(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy) # C + d/dY* B1(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy) # C - F1(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy) # C C21(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy)*d(U)/dT # C + C22(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy)*d(V)/dT # C + C23(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy)*d(W)/dT = # C d/dX* A2(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy) # C + d/dY* B2(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy) # C - F2(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy) # C C31(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy)*d(U)/dT # C + C32(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy)*d(V)/dT # C + C33(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy)*d(W)/dT = # C d/dX* A3(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy) # C + d/dY* B3(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy) # C - F3(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy) # C # C with 'fixed' boundary conditions: # C # C U = FB1(X,Y,T) # C V = FB2(X,Y,T) # C W = FB3(X,Y,T) # C # C or 'free' boundary conditions: # C # C A1*nx + B1*ny = GB1(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy) # C A2*nx + B2*ny = GB2(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy) # C A3*nx + B3*ny = GB3(X,Y,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy) # C # C and initial conditions: # C # C U = U0(X,Y) at T=T0 # C V = V0(X,Y) # C W = W0(X,Y) # C # C where U(X,Y,T),V(X,Y,T),W(X,Y,T) are the unknowns and C11,C12, # C C13,F1,A1,B1,C21,C22,C23,F2,A2,B2,C31,C32,C33,F3,A3,B3,FB1,FB2,FB3, # C GB1,GB2,GB3,U0,V0,W0 are user-supplied functions. # C # C note: # C (nx,ny) = unit outward normal to the boundary # C Ux = d(U)/dX Vx = d(V)/dX Wx = d(W)/dX # C Uy = d(U)/dY Vy = d(V)/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.U F1.Ux F1.Uy F1.V F1.Vx F1.Vy .. +# C + A1.U A1.Ux A1.Uy A1.V A1.Vx A1.Vy .. +# C + B1.U B1.Ux B1.Uy B1.V B1.Vx B1.Vy .. +# C + F2.U F2.Ux F2.Uy F2.V F2.Vx F2.Vy .. +# C + A2.U A2.Ux A2.Uy A2.V A2.Vx A2.Vy .. +# C + B2.U B2.Ux B2.Uy B2.V B2.Vx B2.Vy .. +# C + . . . . . . +# C + . . . . . . +# C + +# C + C11 C12 C13 +# C + C21 C22 C23 +# C + C31 C32 C33 +# C + and +# C + GB1.U GB1.V GB1.W +# C + GB2.U GB2.V GB2.W +# C + GB3.U GB3.V GB3.W +# C + +# C + is always symmetric, where F1.U means d(F1)/d(U), and similarly +# C + for the other terms. In addition, GB1,GB2,GB3 must not depend on +# C + Ux,Uy,Vx,Vy,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 = 4 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 = .TRUE. C GRIDID = .FALSE. IF FINITE ELEMENT GRID CHANGES BETWEEN DUMP, RESTART GRIDID = .FALSE. 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 IMMEDIATELY BELOW: call dtdpx2(nx,ny,xa,xb,ya,yb,hx8z,hy8z,xout8z,yout8z,npts8z) call dtdpzz(ntf,ideg,isolve,symm,neqn,ii8z,ir8z) if (iiwk8z.gt.1) ii8z = iiwk8z if (irwk8z.gt.1) ir8z = irwk8z C *******allocate workspace 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 *******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 if (mod(iprob,4).ne.0) go to 78755 call postpr(tout8z,nsave,xout8z,yout8z,inrg8z,nx,ny,uout,neqn) C *******SURFACE PLOT 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 B1 # C 4 V # C 5 A2 # C 6 B2 # C 7 W # C 8 A3 # C 9 B3 # C############################################################################## IVAR = 1 ivara8z = mod(ivar-1,3)+1 ivarb8z = (ivar-1)/3+1 C############################################################################## C If you don't want to read the FINE PRINT, default ISET1,ISET2,ISINC. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + The tabular output or plots will be made at times: +# C + T(K) = T0 + K*(TF-T0)/NSAVE +# C + for K = ISET1, ISET1+ISINC, ISET1+2*ISINC,..., ISET2 +# C + Enter values for ISET1, ISET2 and ISINC. +# C + +# C + The default is ISET1=0, ISET2=NSAVE, ISINC=1, that is, the tabular +# C + output or plots will be made at all time values for which the +# C + solution has been saved. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## ISET1 = 1 ISET2 = NSAVE ISINC = 1 C############################################################################## C Enter the view latitude, VLAT, and the view longitude, VLON, desired # C for this plot, in degrees. VLAT and VLON must be between 10 and 80 # C degrees; each defaults to 45 degrees. VLAT and VLON are usually # C defaulted. # C############################################################################## VLON = 45.0 VLAT = 45.0 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 UMIN = & alow UMAX = & ahigh 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 ' call dtdprx(tout8z,nsave,iset1,iset2,isinc) do 78756 is8z=iset1,iset2,isinc call dtdpld(xout8z,yout8z,uout(0,0,ivara8z,ivarb8z,is8z),nx,ny,inr &g8z,title,vlon,vlat,umin,umax,tout8z(is8z)) 78756 continue go to 78755 C *******CONTOUR PLOT 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 A # C 3 B # C############################################################################## IVAR = 1 ivara8z = mod(ivar-1,3)+1 ivarb8z = (ivar-1)/3+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 = 'U ' call dtdprx(tout8z,nsave,iset1,iset2,isinc) do 78766 is8z=iset1,iset2,isinc call dtdplg(uout(0,0,ivara8z,ivarb8z,is8z),nx,ny,xa,ya,hx8z,hy8z,i &nrg8z,xbd8z,ybd8z,nbd8z,title,umin,umax,nodist,fillin,tout8z(is8z) &) 78766 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) x = 0.0 y = 0.0 return end subroutine dis8z(x,y,ktri,triden,shape) implicit double precision (a-h,o-z) logical adapt common /parm8z/ pi,r1i,r1r,r2i,r2r,x01,y01,x02,y02 C############################################################################## C 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 TRIDEN may also be a function of the initial triangle number KTRI. # C############################################################################## TRIDEN = 1.0 triden = 1.0/((x-x01)**2+(y-y01)**2) + & 1.0/((x-x02)**2+(y-y02)**2) C############################################################################## C Do you want the triangulation to be graded "adaptively"? # C # C If you answer "yes", 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 two or more times, possibly increasing NTF each time. On the # C first run, the triangulation will be graded as guided by TRIDEN(X,Y), # C but on each subsequent run, information output by the previous run to # C "pde2d.adp" will be used to guide the grading of the new triangulation. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If ADAPT=.TRUE., 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 + 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############################################################################## ADAPT = .TRUE. C EXAG = 0.8 if (adapt) triden = dtdpgr2(x,y,triden**(1.0/exag))**exag 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 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,r1i,r1r,r2i,r2r,x01,y01,x02,y02 U = uu8z(1, 1) Ux = uu8z(2, 1) Uy = uu8z(3, 1) Unorm = Ux*normx + Uy*normy V = uu8z(1, 2) Vx = uu8z(2, 2) Vy = uu8z(3, 2) Vnorm = Vx*normx + Vy*normy W = uu8z(1, 3) Wx = uu8z(2, 3) Wy = uu8z(3, 3) 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,U,Ux,Uy,V,Vx,Vy,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 tr = true(x,y,t,trx,trxx,trxxx) if (kint8z.eq. 1) yd8z = & tr if (kint8z.eq. 2) yd8z = & tr**2 if (kint8z.eq.3) yd8z = U if (kint8z.eq.4) yd8z = U**2 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,U,Ux,Uy,V,Vx,Vy,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 Unorm, # C Vnorm,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,T,U,Ux,Uy,V,Vx,Vy,W,Wx,Wy # C # C They may also be functions of the initial triangle number KTRI. # C # C Recall that the PDEs have the form # C # C C11*d(U)/dT + C12*d(V)/dT + C13*d(W)/dT = d/dX*A1 + d/dY*B1 - F1 # C C21*d(U)/dT + C22*d(V)/dT + C23*d(W)/dT = d/dX*A2 + d/dY*B2 - F2 # C C31*d(U)/dT + C32*d(V)/dT + C33*d(W)/dT = d/dX*A3 + d/dY*B3 - F3 # C # C############################################################################## if (j8z.eq.0) then yd8z = 0.0 C C(1,1) DEFINED if (i8z.eq. -101) yd8z = & 0 C C(1,2) DEFINED if (i8z.eq. -102) yd8z = & 0 C C(1,3) DEFINED if (i8z.eq. -103) yd8z = & 0 C F1 DEFINED if (i8z.eq. 1) yd8z = & V-Ux C A1 DEFINED if (i8z.eq. 2) yd8z = & 0 C B1 DEFINED if (i8z.eq. 3) yd8z = & 0 C C(2,1) DEFINED if (i8z.eq. -201) yd8z = & 0 C C(2,2) DEFINED if (i8z.eq. -202) yd8z = & 0 C C(2,3) DEFINED if (i8z.eq. -203) yd8z = & 0 C F2 DEFINED if (i8z.eq. 4) yd8z = & W C A2 DEFINED if (i8z.eq. 5) yd8z = & Ux C B2 DEFINED if (i8z.eq. 6) yd8z = & 0 C C(3,1) DEFINED if (i8z.eq. -301) yd8z = & 0 C C(3,2) DEFINED if (i8z.eq. -302) yd8z = & 1 C C(3,3) DEFINED if (i8z.eq. -303) yd8z = & 0 C F3 DEFINED if (i8z.eq. 7) yd8z = & 6*V*V+6*U*W C A3 DEFINED if (i8z.eq. 8) yd8z = & -Wx C B3 DEFINED if (i8z.eq. 9) yd8z = & 3*Uy else endif endif return end function u8z(i8z,x,y,ktri,t0) implicit double precision (a-h,o-z) common /parm8z/ pi,r1i,r1r,r2i,r2r,x01,y01,x02,y02 double precision k1I,k1R,k2I,k2R,N1,N1x,N1xx,N2,N2x,N2xx u8z = 0.0 C############################################################################## C Now the initial values must be defined using FORTRAN expressions. # C They may be functions of X and Y, and of the initial triangle number # C KTRI. They may also reference the initial time T0. # C############################################################################## c k1I = r1i c k1R = r1r c k2I = r2i c k2R = r2r c D1 = 4*(x-x01-2*k1I*(y-y01))**2 + 16*k1I**2*(y-y01)**2 + 1/k1I**2 c D1x =8*(x-x01-2*k1I*(y-y01)) c D1xx = 8 c N1 = -4*(x-x01-2*k1I*(y-y01))**2 + 16*k1I**2*(y-y01)**2 + 1/k1I**2 c N1x =-8*(x-x01-2*k1I*(y-y01)) c N1xx = -8 c D2 = 4*(x-x02-2*k2I*(y-y02))**2 + 16*k2I**2*(y-y02)**2 + 1/k2I**2 c D2x =8*(x-x02-2*k2I*(y-y02)) c D2xx = 8 c N2 = -4*(x-x02-2*k2I*(y-y02))**2 + 16*k2I**2*(y-y02)**2 + 1/k2I**2 c N2x =-8*(x-x02-2*k2I*(y-y02)) c N2xx = -8 Tr = true(x,y,0.d0,Trx,Trxx,Trxxx) C U0 DEFINED if (i8z.eq. 1) u8z = & tr c & 16*N1/D1**2 + 16*N2/D2**2 C V0 DEFINED if (i8z.eq. 2) u8z = & trx c & 16*(D1*N1x-2*N1*D1x)/D1**3 + 16*(D2*N2x-2*N2*D2x)/D2**3 C W0 DEFINED if (i8z.eq. 3) u8z = & trxx c & 16*(6*N1*D1x**2+D1**2*N1xx-2*N1*D1*D1xx-4*D1*N1x*D1x)/D1**4 c & + 16*(6*N2*D2x**2+D2**2*N2xx-2*N2*D2*D2xx-4*D2*N2x*D2x)/D2**4 return end function fb8z(i8z,iarc8z,ktri,s,x,y,t) implicit double precision (a-h,o-z) common /parm8z/ pi,r1i,r1r,r2i,r2r,x01,y01,x02,y02 fb8z = 0.0 C NO BOUNDARY CONDITIONS DEFINED ON NEGATIVE ARCS. C TO ADD BCs FOR NEGATIVE ARCS, USE BLOCK BELOW AS MODEL IARC = -1 if (iarc8z.eq.iarc) then C############################################################################## C Enter FORTRAN expressions to define FB1,FB2,FB3 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 U = FB1 # C V = FB2 # C W = FB3 # C # C############################################################################## C reset IXARC(*),IYARC(*) = -1 to use exact values on boundary tr = true(x,y,t,trx,trxx,trxxx) C FB1 DEFINED if (i8z.eq. 1) fb8z = & tr C FB2 DEFINED if (i8z.eq. 2) fb8z = & trx C FB3 DEFINED if (i8z.eq. 3) fb8z = & trxx 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,r1i,r1r,r2i,r2r,x01,y01,x02,y02 zero(f8z) = bign8z*f8z U = uu8z( 1,1) Ux = uu8z( 1,2) Uy = uu8z( 1,3) V = uu8z( 2,1) Vx = uu8z( 2,2) Vy = uu8z( 2,3) W = uu8z( 3,1) Wx = uu8z( 3,2) Wy = uu8z( 3,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,U,Ux,Uy,V,Vx,Vy,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 A3*nx+B3*ny = GB3 # 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 = & 0 C GB3 DEFINED if (i8z.eq. 3) gd8z = & 0 else endif endif C############################################################################## C Enter the arc number (IARC) of a boundary arc with positive arc number. # C############################################################################## IARC = 2 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,U,Ux,Uy,V,Vx,Vy,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 A3*nx+B3*ny = GB3 # 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 = & 0 C GB3 DEFINED if (i8z.eq. 3) gd8z = & 0 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,r1i,r1r,r2i,r2r,x01,y01,x02,y02 U = uu8z(1, 1) Ux = uu8z(2, 1) Uy = uu8z(3, 1) a1 = upr8z(2, 1) b1 = upr8z(3, 1) V = uu8z(1, 2) Vx = uu8z(2, 2) Vy = uu8z(3, 2) a2 = upr8z(2, 2) b2 = upr8z(3, 2) W = uu8z(1, 3) Wx = uu8z(2, 3) Wy = uu8z(3, 3) a3 = upr8z(2, 3) b3 = upr8z(3, 3) 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,B1,V,A2,B2,W,A3,B3 at +# C + the output points. If different variables are to be saved (for +# C + later printing or plotting) the following functions can be used to +# C + re-define the output variables: +# C + define UPRINT(1) to replace U +# C + APRINT(1) A1 +# C + BPRINT(1) B1 +# C + UPRINT(2) V +# C + APRINT(2) A2 +# C + BPRINT(2) B2 +# C + UPRINT(3) W +# C + APRINT(3) A3 +# C + BPRINT(3) B3 +# C + Each function may be a function of +# C + +# C + X,Y,U,Ux,Uy,A1,B1,V,Vx,Vy,A2,B2, +# C + W,Wx,Wy,A3,B3 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(*),BPRINT(*) HERE: c uprint(1) = true(x,y,t,trx,trxx,trxxx) 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,r1i,r1r,r2i,r2r,x01,y01,x02,y02 common /dtdp27/ itask,npes,icomm common /dtdp46/ eps8z,cgtl8z,npmx8z,itype,near8z data lun,lud/0,47/ if (itask.gt.0) return C UOUT(I,J,IDER,IEQ,L) = U_IEQ, if IDER=1 C A_IEQ, if IDER=2 C B_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 COMMENTS TO ACTIVATE) if (lun.eq.0) then lun = 46 open (lun,file='pde2d.m') open (lud,file='pde2d.rdm') write (lun,*) 'fid = fopen(''pde2d.rdm'');' endif do 78753 l=0,nsave if (tout(l).ne.dtdplx(2)) nsave0 = l 78753 continue write (lud,78754) nsave0 write (lud,78754) neqn write (lud,78754) nx write (lud,78754) ny 78754 format (i8) do 78756 i=0,nx do 78755 j=0,ny write (lud,78762) xout(i,j),yout(i,j) 78755 continue 78756 continue do 78761 l=0,nsave0 write (lud,78762) tout(l) do 78760 ieq=1,neqn do 78759 ider=1,3 do 78758 i=0,nx do 78757 j=0,ny if (inrg(i,j).eq.1) then write (lud,78762) uout(i,j,ider,ieq,l) else write (lud,78763) endif 78757 continue 78758 continue 78759 continue 78760 continue 78761 continue 78762 format (e16.8) 78763 format ('NaN') write (lun,*) '% Read solution from pde2d.rdm' write (lun,*) 'NSAVE = fscanf(fid,''%g'',1);' write (lun,*) 'NEQN = fscanf(fid,''%g'',1);' write (lun,*) 'NX = fscanf(fid,''%g'',1);' write (lun,*) 'NY = fscanf(fid,''%g'',1);' write (lun,*) 'if (NX<1 || NY<1)' write (lun,*) ' error(''NX or NY < 1'')' write (lun,*) 'end' if (itype.eq.2) then write (lun,*) 'L0 = 0;' else write (lun,*) 'L0 = 1;' endif write (lun,*) 'T = zeros(NSAVE+1,1);' write (lun,*) 'X = zeros(NY+1,NX+1);' write (lun,*) 'Y = zeros(NY+1,NX+1);' write (lun,*) 'U = zeros(NY+1,NX+1,NSAVE+1,3,NEQN);' write (lun,*) 'for i=0:NX' write (lun,*) 'for j=0:NY' write (lun,*) ' X(j+1,i+1) = fscanf(fid,''%g'',1);' write (lun,*) ' Y(j+1,i+1) = fscanf(fid,''%g'',1);' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'for l=0:NSAVE' write (lun,*) 'T(l+1) = fscanf(fid,''%g'',1);' write (lun,*) 'for ieq=1:NEQN' write (lun,*) 'for ider=1:3' write (lun,*) 'for i=0:NX' write (lun,*) 'for j=0:NY' write (lun,*) & ' U(j+1,i+1,l+1,ider,ieq) = fscanf(fid,''%g'',1);' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) 'xmin = min(min(X(:,:)));' write (lun,*) 'xmax = max(max(X(:,:)));' write (lun,*) 'ymin = min(min(Y(:,:)));' write (lun,*) 'ymax = max(max(Y(:,:)));' write (lun,*) 'hx = 0.1*(xmax-xmin);' write (lun,*) 'hy = 0.1*(ymax-ymin);' write (lun,*) '% Surface plots of each variable' write (lun,*) 'for IEQ=1:NEQN' write (lun,*) '% Plot U_IEQ, if IDER=1' write (lun,*) '% A_IEQ, if IDER=2' write (lun,*) '% B_IEQ, if IDER=3' write (lun,*) 'IDER = 1;' write (lun,*) & 'umin = min(min(min(U(:,:,L0+1:NSAVE+1,IDER,IEQ))));' write (lun,*) & 'umax = max(max(max(U(:,:,L0+1:NSAVE+1,IDER,IEQ))));' write (lun,*) 'if (umax==umin); umax = umin+1; end' write (lun,*) 'for L=L0:NSAVE' write (lun,*) ' figure' write (lun,*) ' surf(X,Y,U(:,:,L+1,IDER,IEQ))' write (lun,*) ' axis([xmin xmax ymin ymax umin umax])' write (lun,*) ' xlabel(''X'')' write (lun,*) ' ylabel(''Y'')' write (lun,*) ' zlabel([''U'',num2str(IEQ)])' write (lun,*) ' title(['' T = '',num2str(T(L+1))])' write (lun,*) ' view(-37.5,30.0)' write (lun,*) 'end' write (lun,*) 'end' write (lun,*) '% ******* Choose variables for vector plots' write (lun,*) '% IEQi,IDERi indicates:' write (lun,*) '% U_IEQi, if IDERi=1' write (lun,*) '% A_IEQi, if IDERi=2' write (lun,*) '% B_IEQi, if IDERi=3' write (lun,*) 'IEQ1 = 1;' write (lun,*) 'IDER1 = 2;' write (lun,*) 'IEQ2 = 1;' write (lun,*) 'IDER2 = 3;' write (lun,*) 'for L=L0:NSAVE' write (lun,*) & ' Umax = max(max(abs(U(:,:,L+1,IDER1,IEQ1))));' write (lun,*) & ' Vmax = max(max(abs(U(:,:,L+1,IDER2,IEQ2))));' write (lun,*) ' figure' write (lun,*) '% For streamlines, replace' write (lun,*) '% quiver with streamslice' write (lun,*) ' quiver(X,Y,U(:,:,L+1,IDER1,IEQ1), ...' write (lun,*) ' U(:,:,L+1,IDER2,IEQ2))' write (lun,*) ' axis([xmin-hx xmax+hx ymin-hy ymax+hy])' write (lun,*) ' xlabel(''X'')' write (lun,*) ' ylabel(''Y'')' write (lun,*) &' title(['' T = '',num2str(T(L+1)),''; Vector Mag. = ('',...' write (lun,*) &' num2str(Umax,''%9.2e''),'','',num2str(Vmax,''%9.2e''),'')'',])' write (lun,*) 'end' return end function true(x,y,t,truex,truexx,truexxx) implicit double precision (a-h,o-z) complex*16 CJ,A(4,4),Ax(4,4),Phi,Phix,Phixx,Phixxx,Phixxxx double precision k1I,k1R,k2I,k2R,kR,kI1,kI2 integer ind1(24),ind2(24),ind3(24),ind4(24),s(24) common /parm8z/ pi,r1i,r1r,r2i,r2r,x01,y01,x02,y02 data ind1/1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4/ data ind2/2,2,3,3,4,4,1,1,3,3,4,4,1,1,2,2,4,4,1,1,2,2,3,3/ data ind3/3,4,2,4,2,3,3,4,1,4,1,3,2,4,1,4,1,2,2,3,1,3,1,2/ data ind4/4,3,4,2,3,2,4,3,4,1,3,1,4,2,4,1,2,1,3,2,3,1,2,1/ data s & /1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1/ C initial conditions for oblique collision k1I = r1i k1R = r1r k2I = r2i k2R = r2r G1I = 2*y01*k1I G2I = 2*y02*k2I G1R = -x01 + k1R/k1I*G1I G2R = -x02 + k2R/k2I*G2I do 10 i=1,4 do 5 j=1,4 Ax(i,j) = 0.0 5 continue Ax(i,i) = 1.0 10 continue CJ = dcmplx(0.d0,1.d0) X1 = x - 2*k1R*y + 12*(k1R**2-k1I**2)*t + G1R Y1 = 2*k1I*(y-12*k1R*t-G1I/(2*k1I)) X2 = x - 2*k2R*y + 12*(k2R**2-k2I**2)*t + G2R Y2 = 2*k2I*(y-12*k2R*t-G2I/(2*k2I)) kR = k1R-k2R kI1 = k1I-k2I kI2 = k1I+k2I A(1,1) = dcmplx(X1,-Y1) A(1,2) = -0.5/k1I A(1,3) = -CJ/dcmplx(kR,kI1) A(1,4) = -CJ/dcmplx(kR,kI2) A(2,1) = 0.5/k1I A(2,2) = dcmplx(X1,Y1) A(2,3) = -CJ/dcmplx(kR,-kI2) A(2,4) = -CJ/dcmplx(kR,-kI1) A(3,1) = CJ/dcmplx(kR,kI1) A(3,2) = CJ/dcmplx(kR,-kI2) A(3,3) = dcmplx(X2,-Y2) A(3,4) = -0.5/k2I A(4,1) = CJ/dcmplx(kR,kI2) A(4,2) = CJ/dcmplx(kR,-kI1) A(4,3) = 0.5/k2I A(4,4) = dcmplx(X2,Y2) Phi = 0.0 Phix = 0.0 Phixx = 0.0 Phixxx = 0.0 Phixxxx = 0.0 do 20 i=1,24 i1 = ind1(i) i2 = ind2(i) i3 = ind3(i) i4 = ind4(i) Phi = Phi + s(i)*A(1,i1)*A(2,i2)*A(3,i3)*A(4,i4) Phix = Phix + s(i)*( & Ax(1,i1)* A(2,i2)* A(3,i3)* A(4,i4) & + A(1,i1)*Ax(2,i2)* A(3,i3)* A(4,i4) & + A(1,i1)* A(2,i2)*Ax(3,i3)* A(4,i4) & + A(1,i1)* A(2,i2)* A(3,i3)*Ax(4,i4)) Phixx = Phixx + 2*s(i)*( & Ax(1,i1)*Ax(2,i2)* A(3,i3)* A(4,i4) & + Ax(1,i1)* A(2,i2)*Ax(3,i3)* A(4,i4) & + Ax(1,i1)* A(2,i2)* A(3,i3)*Ax(4,i4) & + A(1,i1)*Ax(2,i2)*Ax(3,i3)* A(4,i4) & + A(1,i1)*Ax(2,i2)* A(3,i3)*Ax(4,i4) & + A(1,i1)* A(2,i2)*Ax(3,i3)*Ax(4,i4)) Phixxx = Phixxx + 6*s(i)*( & A(1,i1)*Ax(2,i2)*Ax(3,i3)*Ax(4,i4) & + Ax(1,i1)* A(2,i2)*Ax(3,i3)*Ax(4,i4) & + Ax(1,i1)*Ax(2,i2)* A(3,i3)*Ax(4,i4) & + Ax(1,i1)*Ax(2,i2)*Ax(3,i3)* A(4,i4)) Phixxxx = Phixxxx + 24*s(i)*Ax(1,i1)*Ax(2,i2)*Ax(3,i3)*Ax(4,i4) 20 continue true = 2*dreal(Phixx/Phi - (Phix/Phi)**2) truex= 2*dreal(Phixxx/Phi - 3*Phix*Phixx/Phi**2 + 2*(Phix/Phi)**3) truexx = 2*dreal(Phixxxx/Phi - 4*Phix*Phixxx/Phi**2 & - 3*(Phixx/Phi)**2 + 12*Phix**2*Phixx/Phi**3 - 6*(Phix/Phi)**4) c truexxx = 2*dreal(-2*Phix*Phixxxx/Phi**2 - 7*Phixx*Phixxx/Phi**2 c & + 14*Phix**2*Phixxx/Phi**3 + 30*Phix*Phixx**2/Phi**3 c & - 60*Phix**3*Phixx/Phi**4 + 24*(Phix/Phi)**5) return end