C **************************
C * PDE2D 9.5 MAIN PROGRAM *
C **************************
C##############################################################################
C##############################################################################
C *** 2D PROBLEM SOLVED (GALERKIN METHOD) ***
C##############################################################################
C Is double precision mode to be used? Double precision is recommended. #
C #
C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C + If double precision mode is used, variables and functions assigned +#
C + names beginning with a letter in the range A-H or O-Z will be DOUBLE +#
C + PRECISION, and you should use double precision constants and FORTRAN +#
C + expressions throughout; otherwise such variables and functions will +#
C + be of type REAL. In either case, variables and functions assigned +#
C + names beginning with I,J,K,L,M or N will be of INTEGER type. +#
C + +#
C + It is possible to convert a single precision PDE2D program to double +#
C + precision after it has been created, using an editor. Just change +#
C + all occurrences of "real" to "double precision" +#
C + " tdp" to "dtdp" (note leading blank) +#
C + Any user-written code or routines must be converted "by hand", of +#
C + course. To convert from double to single, reverse the changes. +#
C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
implicit double precision (a-h,o-z)
parameter (neqnmx= 99)
parameter (ndelmx= 20)
parameter (nv0=0,nt0=0,npgrid=0,nqgrid=0)
C##############################################################################
C NXGRID = number of X-grid lines #
C##############################################################################
PARAMETER (NXGRID = 4)
C##############################################################################
C NYGRID = number of Y-grid lines #
C##############################################################################
PARAMETER (NYGRID = 4)
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= 5000000)
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 + +#
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 = 10)
PARAMETER (NY = 10)
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 = 380)
common/parm8z/ pi,W ,AMU
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 = 1
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 + 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 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
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 ++++++++++++++++++++++++++ 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##############################################################################
W =
& 2*PI*3.5269 D14
AMU =
& 12.566 D-7
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) =
& -4.D-6
XGRID(2) =
& -1.D-6
XGRID(3) =
& 1.D-6
XGRID(NXGRID) =
& 4.D-6
C YGRID DEFINED
YGRID(1) =
& -4.D-6
YGRID(2) =
& -1.D-6
YGRID(3) =
& 1.D-6
YGRID(NYGRID) =
& 4.D-6
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) =
& 1
IXARC(2) =
& 1
IYARC(1) =
& 1
IYARC(2) =
& 1
C##############################################################################
C How many triangles (NTF) are desired for the final triangulation? #
C##############################################################################
NTF = 5000
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 = 2
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 =
& 1.13 D7
TF =
& 1.14 D7
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 = .FALSE.
TOLER(1) = 0.01
C##############################################################################
C If you don't want to read the FINE PRINT, it is safe (though possibly #
C very inefficient) to enter 'no'. #
C #
C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C + If your time-dependent problem is linear with all PDE and boundary +#
C + condition coefficients independent of time except inhomogeneous +#
C + terms, then a large savings in execution time may be possible if +#
C + this is recognized (the LU decomposition computed on the first step +#
C + can be used on subsequent steps). Is this the case for your +#
C + problem? (Caution: if you answer 'yes' when you should not, you +#
C + will get incorrect results with no warning.) +#
C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
NOUPDT = .FALSE.
C##############################################################################
C The time stepsize will be constant, DT = (TF-T0)/NSTEPS. Enter a #
C value for NSTEPS, the number of time steps. #
C #
C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C + If you later turn on adaptive step control, the time stepsize will be+#
C + chosen adaptively, between an upper limit of DTMAX = (TF-T0)/NSTEPS +#
C + and a lower limit of 0.0001*DTMAX. +#
C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
NSTEPS =
& NSAVE
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 = .FALSE.
C##############################################################################
C PDE2D solves the system of equations: #
C #
C C11(X,Y,T,H,Hx,Hy,E,Ex,Ey)*d(H)/dT #
C + C12(X,Y,T,H,Hx,Hy,E,Ex,Ey)*d(E)/dT = #
C d/dX* A1(X,Y,T,H,Hx,Hy,E,Ex,Ey) #
C + d/dY* B1(X,Y,T,H,Hx,Hy,E,Ex,Ey) #
C - F1(X,Y,T,H,Hx,Hy,E,Ex,Ey) #
C C21(X,Y,T,H,Hx,Hy,E,Ex,Ey)*d(H)/dT #
C + C22(X,Y,T,H,Hx,Hy,E,Ex,Ey)*d(E)/dT = #
C d/dX* A2(X,Y,T,H,Hx,Hy,E,Ex,Ey) #
C + d/dY* B2(X,Y,T,H,Hx,Hy,E,Ex,Ey) #
C - F2(X,Y,T,H,Hx,Hy,E,Ex,Ey) #
C #
C with 'fixed' boundary conditions: #
C #
C H = FB1(X,Y,T) #
C E = FB2(X,Y,T) #
C #
C or 'free' boundary conditions: #
C #
C A1*nx + B1*ny = GB1(X,Y,T,H,Hx,Hy,E,Ex,Ey) #
C A2*nx + B2*ny = GB2(X,Y,T,H,Hx,Hy,E,Ex,Ey) #
C #
C and initial conditions: #
C #
C H = H0(X,Y) at T=T0 #
C E = E0(X,Y) #
C #
C where H(X,Y,T) and E(X,Y,T) are the unknowns and C11,C12,F1,A1,B1, #
C C21,C22,F2,A2,B2,FB1,FB2,GB1,GB2,H0,E0 are user-supplied functions. #
C #
C note: #
C (nx,ny) = unit outward normal to the boundary #
C Hx = d(H)/dX Ex = d(E)/dX #
C Hy = d(H)/dY Ey = d(E)/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.H F1.Hx F1.Hy F1.E F1.Ex F1.Ey +#
C + A1.H A1.Hx A1.Hy A1.E A1.Ex A1.Ey +#
C + B1.H B1.Hx B1.Hy B1.E B1.Ex B1.Ey +#
C + F2.H F2.Hx F2.Hy F2.E F2.Ex F2.Ey +#
C + A2.H A2.Hx A2.Hy A2.E A2.Ex A2.Ey +#
C + B2.H B2.Hx B2.Hy B2.E B2.Ex B2.Ey +#
C + +#
C + C11 C12 +#
C + C21 C22 +#
C + and +#
C + GB1.H GB1.E +#
C + GB2.H GB2.E +#
C + +#
C + is always symmetric, where F1.H means d(F1)/d(H), and similarly +#
C + for the other terms. In addition, GB1,GB2 must not depend on +#
C + Hx,Hy,Ex,Ey. +#
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.
FDIFF = .FALSE.
C##############################################################################
C You may calculate one or more integrals (over the entire region) of #
C some functions of the solution and its derivatives. How many integrals #
C (NINT), if any, do you want to calculate? #
C #
C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C + In the FORTRAN program created by the preprocessor, the computed +#
C + values of the integrals will be returned in the vector SINT8Z. If +#
C + several iterations or time steps are done, only the last computed +#
C + values are saved in SINT8Z (all values are printed). +#
C + +#
C + A limiting value, SLIM8Z(I), for the I-th integral can be set +#
C + below in the main program. The computations will then stop +#
C + gracefully whenever SINT8Z(I) > SLIM8Z(I), for any I=1...NINT. +#
C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
NINT = 1
C##############################################################################
C You may calculate one or more boundary integrals (over the entire #
C boundary) of some functions of the solution and its derivatives. How #
C many boundary integrals (NBINT), if any, do you want to calculate? #
C #
C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C + In the FORTRAN program created by the preprocessor, the computed +#
C + values of the integrals will be returned in the vector BINT8Z. If +#
C + several iterations or time steps are done, only the last computed +#
C + values are saved in BINT8Z (all values are printed). +#
C + +#
C + A limiting value, BLIM8Z(I), for the I-th boundary integral can be +#
C + set below in the main program. The computations will then stop +#
C + gracefully whenever BINT8Z(I) > BLIM8Z(I), for any I=1...NBINT. +#
C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
NBINT = 0
ndel = 0
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 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 = .FALSE.
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
call postpr(tout8z,nsave,xout8z,yout8z,inrg8z,nx,ny,uout,neqn)
C *******1D CROSS-SECTION PLOTS
C##############################################################################
C Enter a value for IVAR, to select the variable to be plotted or #
C printed: #
C IVAR = 1 means H (possibly as modified by UPRINT,..) #
C 2 A1 #
C 3 B1 #
C 4 E #
C 5 A2 #
C 6 B2 #
C##############################################################################
IVAR = 2
C T IS VARIABLE
ics8z = 3
C##############################################################################
C One-dimensional plots of the output variable as a function of T will #
C be made, at the output grid points (X,Y) closest to #
C (XCROSS(I),YCROSS(J)), I=1,...,NXVALS, J=1,...,NYVALS #
C #
C Enter values for NXVALS, XCROSS(1),...,XCROSS(NXVALS), #
C and NYVALS, YCROSS(1),...,YCROSS(NYVALS) #
C##############################################################################
NXVALS = 1
XCROSS(1) =
& 0
NYVALS = 1
YCROSS(1) =
& 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
C##############################################################################
C Enter a title, WITHOUT quotation marks. A maximum of 40 characters #
C are allowed. The default is no title. #
C##############################################################################
TITLE = ' '
TITLE = 'Solution norm vs. beta '
is8z = 0
do 78758 ixv8z=1,nxvals
do 78757 jyv8z=1,nyvals
call dtdpzx(xcross(ixv8z),xa,xb,nx,ix8z)
call dtdpzx(ycross(jyv8z),ya,yb,ny,jy8z)
call dtdplp(ics8z,ivar,tout8z,nsave,xout8z,yout8z,inrg8z,nx,ny,uou
&t,neqn,title,umin,umax,ix8z,jy8z,is8z)
78757 continue
78758 continue
go to 78755
C *******SURFACE PLOT
C##############################################################################
C Enter a value for IVAR, to select the variable to be plotted or #
C printed: #
C IVAR = 1 means H (possibly as modified by UPRINT,..) #
C 2 A1 #
C 3 B1 #
C 4 E #
C 5 A2 #
C 6 B2 #
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 = 0
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
C##############################################################################
C Enter a title, WITHOUT quotation marks. A maximum of 40 characters #
C are allowed. The default is no title. #
C##############################################################################
TITLE = ' '
TITLE = 'H '
call dtdprx(tout8z,nsave,iset1,iset2,isinc)
do 78759 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))
78759 continue
C *******SURFACE PLOT
C##############################################################################
C Enter a value for IVAR, to select the variable to be plotted or #
C printed: #
C IVAR = 1 means H (possibly as modified by UPRINT,..) #
C 2 A1 #
C 3 B1 #
C 4 E #
C 5 A2 #
C 6 B2 #
C##############################################################################
IVAR = 4
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 = 0
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
C##############################################################################
C Enter a title, WITHOUT quotation marks. A maximum of 40 characters #
C are allowed. The default is no title. #
C##############################################################################
TITLE = ' '
TITLE = 'E '
call dtdprx(tout8z,nsave,iset1,iset2,isinc)
do 78760 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))
78760 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 H (possibly as modified by UPRINT,..) #
C 2 A1 #
C 3 B1 #
C 4 E #
C 5 A2 #
C 6 B2 #
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 = 0
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 = .FALSE.
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 = 'H '
call dtdprx(tout8z,nsave,iset1,iset2,isinc)
do 78761 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)
&)
78761 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 H (possibly as modified by UPRINT,..) #
C 2 A1 #
C 3 B1 #
C 4 E #
C 5 A2 #
C 6 B2 #
C##############################################################################
IVAR = 4
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 = 0
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 = .FALSE.
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 = 'E '
call dtdprx(tout8z,nsave,iset1,iset2,isinc)
do 78762 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)
&)
78762 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,W ,AMU
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 =
& EXP(-(X/1.D-6)**2-(Y/1.D-6)**2)
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 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 +#
C + the first run, the triangulation will be graded as guided by +#
C + TRIDEN(X,Y), but on each subsequent run, information output by the +#
C + previous run to "pde2d.adp" will be used to guide the grading of the +#
C + new triangulation. +#
C +#
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. If RESTRT=.TRUE., GRIDID=.FALSE., and +#
C + T0,TF are incremented each pass through the DO loop, it is possible +#
C + in this way to solve a time-dependent problem with an adaptive, +#
C + moving, grid. +#
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 = .FALSE.
C
EXAG = 1.5
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,W ,AMU
H = uu8z(1, 1)
Hx = uu8z(2, 1)
Hy = uu8z(3, 1)
Hnorm = Hx*normx + Hy*normy
E = uu8z(1, 2)
Ex = uu8z(2, 2)
Ey = uu8z(3, 2)
Enorm = Ex*normx + Ey*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,H,Hx,Hy,E,Ex,Ey 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 =
& abs(H)+abs(E)
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,H,Hx,Hy,E,Ex,Ey 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 Hnorm, #
C Enorm. #
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,H,Hx,Hy,E,Ex,Ey #
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(H)/dT + C12*d(E)/dT = d/dX*A1 + d/dY*B1 - F1 #
C C21*d(H)/dT + C22*d(E)/dT = d/dX*A2 + d/dY*B2 - F2 #
C #
C##############################################################################
IF (ABS(X).LE.1.E-6 .AND. ABS(Y).LE.1.E-6) THEN
C GLASS REGION
EPS = 8.854 D-12*1.55**2
ELSE
C SILICA REGION
EPS = 8.854 D-12*1.50**2
ENDIF
BETA = T
DENOM = W**2*AMU*EPS - BETA**2
G1 = SIN(X/0.1E-6 + Y/0.2E-6)
G2 = SIN(X/0.2E-6 + Y/0.1E-6)
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 F1 DEFINED
if (i8z.eq. 1) yd8z =
& -W*AMU*H + G1
C A1 DEFINED
if (i8z.eq. 2) yd8z =
& (W*AMU*Hx - BETA*Ey)/DENOM
C B1 DEFINED
if (i8z.eq. 3) yd8z =
& (W*AMU*Hy + BETA*Ex)/DENOM
C C(2,1) DEFINED
if (i8z.eq. -201) yd8z =
& 0
C C(2,2) DEFINED
if (i8z.eq. -202) yd8z =
& 0
C F2 DEFINED
if (i8z.eq. 4) yd8z =
& -W*EPS*E + G2
C A2 DEFINED
if (i8z.eq. 5) yd8z =
& (W*EPS*Ex + BETA*Hy)/DENOM
C B2 DEFINED
if (i8z.eq. 6) yd8z =
& (W*EPS*Ey - BETA*Hx)/DENOM
else
endif
endif
return
end
function u8z(i8z,x,y,ktri,t0)
implicit double precision (a-h,o-z)
common/parm8z/ pi,W ,AMU
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 H0 DEFINED
C if (i8z.eq. 1) u8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
C E0 DEFINED
C if (i8z.eq. 2) u8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
return
end
function fb8z(i8z,iarc8z,ktri,s,x,y,t)
implicit double precision (a-h,o-z)
common/parm8z/ pi,W ,AMU
fb8z = 0.0
C NO BOUNDARY CONDITIONS DEFINED ON NEGATIVE ARCS.
C TO ADD BCs FOR NEGATIVE ARCS, USE BLOCK BELOW AS MODEL
IARC = 0
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 H = FB1 #
C E = FB2 #
C #
C##############################################################################
C FB1 DEFINED
C if (i8z.eq. 1) fb8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
C FB2 DEFINED
C if (i8z.eq. 2) fb8z =
C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]
endif
return
end
subroutine gb8z(gd8z,i8z,j8z,iarc8z,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,W ,AMU
zero(f8z) = bign8z*f8z
H = uu8z( 1,1)
Hx = uu8z( 1,2)
Hy = uu8z( 1,3)
E = uu8z( 2,1)
Ex = uu8z( 2,2)
Ey = 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,H,Hx,Hy,E,Ex,Ey 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(E)
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,W ,AMU
H = uu8z(1, 1)
Hx = uu8z(2, 1)
Hy = uu8z(3, 1)
a1 = upr8z(2, 1)
b1 = upr8z(3, 1)
E = uu8z(1, 2)
Ex = uu8z(2, 2)
Ey = 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 H,A1,B1,E,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 H +#
C + APRINT(1) A1 +#
C + BPRINT(1) B1 +#
C + UPRINT(2) E +#
C + APRINT(2) A2 +#
C + BPRINT(2) B2 +#
C + Each function may be a function of +#
C + +#
C + X,Y,H,Hx,Hy,A1,B1,E,Ex,Ey,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 H. 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) =
& SINT(1)
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,W ,AMU
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 C! COMMENTS TO ACTIVATE)
C! if (lun.eq.0) then
C! lun = 46
C! open (lun,file='pde2d.m')
C! open (lud,file='pde2d.rdm')
C! write (lun,*) 'fid = fopen(''pde2d.rdm'');'
C! endif
C! do 78753 l=0,nsave
C! if (tout(l).ne.dtdplx(2)) nsave0 = l
C!78753 continue
C! write (lud,78754) nsave0
C! write (lud,78754) neqn
C! 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 pde2d.m
C! call mtdp2dg(itype,lun)
return
end