C ************************** C * PDE2D 9.5 MAIN PROGRAM * C ************************** C *** 1D PROBLEM SOLVED (COLLOCATION 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) dimension saveE(21,5),errE(21,5),rkE(21) parameter (neqnmx= 99) C############################################################################## C NXGRID = number of X-grid lines # C############################################################################## PARAMETER (NXGRID = 101) C############################################################################## C How many differential equations (NEQN) are there in your problem? # C############################################################################## PARAMETER (NEQN = 2) C DIMENSIONS OF WORK ARRAYS C SET TO 1 FOR AUTOMATIC ALLOCATION PARAMETER (IRWK8Z= 1) PARAMETER (IIWK8Z= 1) PARAMETER (NXP8Z=1001,KDEG8Z=1) C############################################################################## C The solution is saved on a uniform grid of NX+1 points # C XA + I*(XB-XA)/NX # C I=0,...,NX. Enter a value for NX (suggested value = 50). # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you want to save the solution at an arbitrary user-specified set +# C + of NX+1 points, enter -NX. +# 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 + collocation or integration point. For most problems this obviously +# C + will produce less accurate output or plots, but for certain (rare) +# C + problems, a solution component may be much less noisy when plotted +# C + only at collocation or integration points. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## PARAMETER (NX = 50) PARAMETER (NSAVE = 1) common/parm8z/ pi,p,rk,a ,b dimension xgrid(nxgrid),xout8z(0:nx),xcross(100),tout8z(0:nsave) C dimension xres8z(nxp8z),ures8z(neqn,nxp8z) allocatable iwrk8z(:),rwrk8z(:) C dimension iwrk8z(iiwk8z),rwrk8z(irwk8z) character*40 title logical linear,crankn,noupdt,nodist,fillin,evcmpx,adapt,plot,lsqfi &t,fdiff,econ8z,ncon8z,restrt,gridid common/dtdp14/ sint8z(20),bint8z(20),slim8z(20),blim8z(20) common/dtdp15/ evlr8z,ev0r,evli8z,ev0i,evcmpx common/dtdp16/ p8z,evr8z(50),evi8z(50) common/dtdp19/ toler(neqnmx),adapt common/dtdp30/ econ8z,ncon8z common/dtdp42/ nxa8z,kd8z common/dtdp43/ work8z(nxp8z+3) common/dtdp45/ perdc(neqnmx) common/dtdp46/ eps8z,cgtl8z,npmx8z,itype,near8z common/dtdp62/ amin8z(2*neqnmx),amax8z(2*neqnmx) common/dtdp75/ nx18z,xa,xb,uout(0:nx,2,neqn,0:nsave) pi = 4.0*atan(1.d0) nxa8z = nxp8z nx18z = nx+1 kd8z = kdeg8z C############################################################################## C If you don't want to read the FINE PRINT, default NPROB. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you want to solve several similar problems in the same run, set +# C + NPROB equal to the number of problems you want to solve. Then NPROB +# C + loops through the main program will be done, with IPROB=1,...,NPROB, +# C + and you can make the problem parameters vary with IPROB. NPROB +# C + defaults to 1. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## NPROB = 21 do 78755 iprob=1,nprob C############################################################################## C PDE2D solves the time-dependent system (note: U,F,G,U0 may be vectors, # C C,RHO may be matrices): # C # C C(X,T,U,Ux)*d(U)/dT = F(X,T,U,Ux,Uxx) # C # C or the steady-state system: # C # C F(X,U,Ux,Uxx) = 0 # C # C or the linear and homogeneous eigenvalue system: # C # C F(X,U,Ux,Uxx) = lambda*RHO(X)*U # C # C with boundary conditions: # C # C G(X,[T],U,Ux) = 0 # C (periodic boundary conditions are also permitted) # C # C at two X values. # C # C For time-dependent problems there are also initial conditions: # C # C U = U0(X) at T=T0 # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If your PDEs involve the solution at points other than X, the +# C + function +# C + (D)OLDSOL1(IDER,IEQ,XX,KDEG) +# C + will interpolate (using interpolation of degree KDEG=1,2 or 3) to XX +# C + the function saved in UOUT(*,IDER,IEQ,ISET) on the last time step or +# C + iteration (ISET) for which it has been saved. Thus, for example, if +# C + IDER=1, this will return the latest value of component IEQ of the +# C + solution at XX, assuming this has not been modified using UPRINT... +# C + If your equations involve integrals of the solution, for example, +# C + you can use (D)OLDSOL1 to approximate these using the solution from +# C + 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 - 1/(1+U**10) = 0, where U = UR + UI*I +# C + would be difficult to split up analytically, but using FORTRAN +# C + expressions it is easy: +# C + F1 = -UIxx - REAL(1.0/(1.0+CMPLX(UR,UI)**10)) +# C + F2 = URxx - 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############################################################################## a = & 1. b = & 0.01 P = 1.5 rk = 2*pi*float(iprob-1)/float(nprob-1)/a C############################################################################## C A collocation finite element method is used, with cubic Hermite # C basis functions on the subintervals defined by the grid points: # C XGRID(1),XGRID(2),...,XGRID(NXGRID) # C You will first be prompted for NXGRID, the number of X-grid points, # C then for XGRID(1),...,XGRID(NXGRID). Any points defaulted will be # C uniformly spaced between the points you define; the first and last # C points cannot be defaulted. The interval over which the PDE system # C is to be solved is then: # C XGRID(1) < X < XGRID(NXGRID) # C # C############################################################################## call dtdpwx(xgrid,nxgrid,0) C XGRID DEFINED XGRID(1) = & 0 XGRID(NXGRID/2) = a-2*b XGRID(NXGRID/2+10) = a-b XGRID(NXGRID) = & a call dtdpwx(xgrid,nxgrid,1) C *******EIGENVALUE PROBLEM C############################################################################## C The shifted inverse power method will be used to find one eigenvalue # C and the corresponding eigenfunction. # C # C You can later find ALL eigenvalues (without eigenfunctions), including # C complex eigenvalues, by changing ITYPE from 3 to 4 below in the FORTRAN # C program. All eigenvalues will be written to a file 'pde2d.eig'; the 50 # C eigenvalues closest to P8Z (P8Z is 0 by default but can be set by the # C user below) will be printed and also available in the main program as # C EVR8Z(k) + i*EVI8Z(k), for k=1,...,min(50,N). # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Caution: for ITYPE=4, the memory and computer time requirements are +# C + much greater than for ITYPE=3, for a given grid. However, this +# C + calculation runs efficiently on multi-processor machines, under MPI. +# C + +# C + The calculation is also dramatically faster for symmetric 1D +# C + Galerkin problems, and symmetric 2D Galerkin problems (IDEG=-1,-3) +# C + because the algorithm can take advantage of the band structure of A +# C + in finding all eigenvalues of the generalized eigenvalue problem +# C + Az=lambda*Bz, but only when A is symmetric and B is diagonal (hence +# C + the requirement that IDEG=-1 or -3 for 2D Galerkin problems; the RHO +# C + matrix must also be diagonal (if NEQN > 1) and of constant sign). +# C + High precision is STRONGLY recommended in these cases. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## ITYPE = 4 C C P8Z = 0.0 linear = .true. t0 = 0.0 dt = 1.0 crankn = .false. noupdt = .true. C############################################################################## C If you don't want to read the FINE PRINT, enter 'no'. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Are your eigenvalue PDEs complex? If 'yes', you must write each of +# C + the unknown eigenfunctions, but NOT the eigenvalue (i.e., pretend +# C + the eigenvalue is real even though it may not be), in terms of its +# C + real and imaginary parts, and separate the equations into real and +# C + imaginary parts. Then it is assumed that your unknowns represent +# C + alternately the real and imaginary parts of the eigenfunctions, and +# C + that your equations represent alternately the real and imaginary +# C + parts of the eigenvalue differential equations. Note that your +# C + first equation must represent exactly the real part of the first +# C + complex PDE, your second equation must represent exactly the +# C + imaginary part, etc. You must not reorder the equations, nor +# C + multiply any equation through by a constant. +# C + +# C + Even if your eigenvalue PDEs are real, if a desired eigenvalue/ +# C + eigenfunction is complex, you should answer 'yes' and treat the +# C + PDEs as if they were complex, as outlined above. If you are going +# C + to set ITYPE=4 to find all eigenvalues without eigenfunctions, +# C + however, you should answer 'no' even if some eigenvalues are complex.+# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C # C If you don't want to read the FINE PRINT above, enter 'no'. # C # C############################################################################## EVCMPX = .FALSE. C############################################################################## C The (shifted) inverse power method will be used to find the eigenvalue # C closest to EV0R; enter a value for EV0R. The closer you choose EV0R # C to the desired eigenvalue, the faster the convergence will be. The # C default is EV0R = 0.0, that is, the smallest eigenvalue (in absolute # C value) is found. # C############################################################################## EV0R = 0.0 EV0R = & 0 C############################################################################## C How many iterations (NSTEPS) of the inverse power method do you want # C to do? NSTEPS defaults to 25. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + After convergence of the inverse power method, the solution each +# C + iteration will be a normalized eigenfunction corresponding to the +# C + eigenvalue closest to EV0R (or EV0R+EV0I*I) and this eigenvalue will +# C + be printed. In the FORTRAN program created by the preprocessor, +# C + the last computed estimate of this eigenvalue will be returned as +# C + EVLR8Z (or EVLR8Z+EVLI8Z*I), and ECON8Z will be .TRUE. if the +# C + inverse power method has converged. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## NSTEPS = 25 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 = 0 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 lsqfit = .false. C############################################################################## C If you don't want to read the FINE PRINT, enter 'no'. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Do you want to read the initial conditions from the restart file, +# C + if it exists (and use the conditions supplied above if it does not +# C + exist)? +# C + +# C + If so, PDE2D will dump the final solution at the end of each run +# C + into a restart file "pde2d.res". Thus the usual procedure for +# C + using this dump/restart option is to make sure there is no restart +# C + file in your directory left over from a previous job, then the +# C + first time you run this job, the initial conditions supplied above +# C + will be used, but on the second and subsequent runs the restart file +# C + from the previous run will be used to define the initial conditions. +# C + +# C + You can do all the "runs" in one program, by setting NPROB > 1. +# C + Each pass through the DO loop, T0,TF,NSTEPS and possibly other +# C + parameters may be varied, by making them functions of IPROB. +# C + +# C + If the 2D or 3D collocation method is used, the coordinate +# C + transformation should not change between dump and restart. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## RESTRT = .FALSE. C GRIDID = .FALSE. IF FINITE ELEMENT GRID CHANGES BETWEEN DUMP, RESTART GRIDID = .TRUE. C############################################################################## C If you do not have periodic boundary conditions, enter IPERDC=0. # C # C Enter IPERDC=1 for periodic conditions at X = XGRID(1),XGRID(NXGRID) # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + When periodic boundary conditions are selected, they apply to all +# C + variables by default. To turn off periodic boundary conditions on +# C + the I-th variable, set PERDC(I) to 0 below in the main program and +# C + set the desired boundary conditions in subroutine GB8Z, "by hand". +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## IPERDC = 1 C############################################################################## C The solution is saved on a uniform grid of NX+1 points, covering the # C interval (XA,XB). Enter values for XA,XB. These variables are usually # C defaulted. # C # C The defaults are XA = XGRID(1), XB = XGRID(NXGRID) # C # C############################################################################## C defaults for xa,xb xa = xgrid(1) xb = xgrid(nxgrid) C DEFINE XA,XB IMMEDIATELY BELOW: call dtdpx1(nx,xa,xb,hx8z,xout8z,npts8z) C SOLUTION SAVED EVERY NOUT ITERATIONS NOUT = NSTEPS call dtdp1q(nxgrid,neqn,ii8z,ir8z) if (iiwk8z.gt.1) ii8z = iiwk8z if (irwk8z.gt.1) ir8z = irwk8z C *******allocate workspace allocate (iwrk8z(ii8z),rwrk8z(ir8z)) C *******DRAW GRID POINTS? PLOT = .FALSE. C *******call pde solver call dtdp1x(xgrid, nxgrid, neqn, nint, nbint, xout8z, uout, tout8z &, iperdc, plot, lsqfit, fdiff, npts8z, t0, dt, nsteps, nout, nsave &, crankn, noupdt, itype, linear, rwrk8z, ir8z, iwrk8z, ii8z, restr &t, gridid) deallocate (iwrk8z,rwrk8z) rkE(iprob) = rk*a do k=1,5 saveE(iprob,k) = EVR8Z(k) al = sqrt(EVR8Z(k)) errE(iprob,k) = cos(al*a) + P*sin(al*a)/(al*a) - cos(rk*a) enddo if (itype.eq.4) go to 78755 C *******UNCOMMENT TO POSTPROCESS ONLY IF CONVERGED C IF (.NOT.ECON8Z) GO TO 78755 C *******read from restart file to array ures8z C call dtdpr1(1,xres8z,nxp8z,ures8z,neqn) C *******write array ures8z back to restart file C call dtdpr1(2,xres8z,nxp8z,ures8z,neqn) C *******call user-written postprocessor call postpr(tout8z,nsave,xout8z,nx,uout,neqn) C *******LINE PLOTS C############################################################################## C Enter a value for IVAR, to select the variable to be plotted or # C printed: # C IVAR = 1 means UR (possibly as modified by UPRINT,..) # C 2 URx # C 3 UI # C 4 UIx # C############################################################################## IVAR = 1 C X IS VARIABLE ics8z = 1 ISET1 = 1 ISET2 = NSAVE ISINC = 1 C alow = amin8z(ivar) ahigh = amax8z(ivar) C############################################################################## C Specify the range (UMIN,UMAX) for the dependent variable axis. UMIN # C and UMAX are often defaulted. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + By default, each plot will be scaled to just fit in the plot area. +# C + For a common scaling, you may want to set UMIN=ALOW, UMAX=AHIGH. +# C + ALOW and AHIGH are the minimum and maximum values over all output +# C + points and over all saved time steps or iterations. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## UMIN = 0.0 UMAX = 0.0 C############################################################################## C Enter a title, WITHOUT quotation marks. A maximum of 40 characters # C are allowed. The default is no title. # C############################################################################## TITLE = ' ' TITLE = 'UR ' call dtdprx(tout8z,nsave,iset1,iset2,isinc) do 78756 is8z=iset1,iset2,isinc call dtdpzp(ics8z,ivar,tout8z,nsave,xout8z,nx,uout,neqn,title,umin &,umax,ix8z,is8z) 78756 continue 78755 continue open (11,file='Kronig.txt') do i=1,nprob write(11,9) rkE(i),(saveE(i,j),j=1,5) 9 format (6F15.5) enddo do i=1,nprob print 9, (errE(i,j),j=1,5) enddo call endgks stop end subroutine pdes8z(yd8z,i8z,j8z,kint8z,x,t,uu8z) implicit double precision (a-h,o-z) parameter (neqnmx= 99) parameter (NEQN= 2) C un8z(1,I),un8z(2,I),... hold the (rarely used) values C of UI,UIx,... from the previous iteration or time step common /dtdp4x/un8z(3,neqnmx) common /dtdp11/normx double precision normx,uu8z(3,neqnmx) dimension rho(neqn,neqn) common/parm8z/ pi,p,rk,a ,b UR = uu8z(1, 1) URx = uu8z(2, 1) URxx= uu8z(3, 1) UI = uu8z(1, 2) UIx = uu8z(2, 2) UIxx= uu8z(3, 2) if (i8z.eq.0) then yd8z = 0.0 C############################################################################## C Enter FORTRAN expressions for the functions whose integrals are to be # C calculated and printed. They may be functions of # C # C X,UR,URx,URxx,UI,UIx,UIxx and (if applicable) T # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you only want to integrate a function over part of the interval, +# C + define that function to be zero on the rest of the interval. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C INTEGRAL1 DEFINED C if (kint8z.eq.1) yd8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] C############################################################################## C Enter FORTRAN expressions for the functions whose "integrals" (sum # C over two boundary points) are to be calculated and printed. They may # C be functions of # C # C X,UR,URx,URxx,UI,UIx,UIxx and (if applicable) T # C # C The unit outward normal, NORMx (=1 at right endpoint, -1 at left), # C may also be referenced. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If you only want to "integrate" a function over one boundary point, +# C + define that function to be zero at the other point. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C BND. INTEGRAL1 DEFINED C if (kint8z.eq.-1) yd8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] else C############################################################################## C Now enter FORTRAN expressions to define the PDE coefficients. # C RHO11,RHO12,RHO21,RHO22 may be functions of X, while F1,F2 # C may be functions of # C # C X,UR,URx,URxx,UI,UIx,UIxx # C # C Recall that the PDEs have the form # C # C F1 = lambda*(RHO11*UR + RHO12*UI) # C F2 = lambda*(RHO21*UR + RHO22*UI) # C # C############################################################################## W = 0 if (x.gt.a-b) W = 2*P/(b*a) if (j8z.eq.0) then yd8z = 0.0 C F1 DEFINED if (i8z.eq. 1) yd8z = & URxx - 2*rk*UIx - rk**2*UR - W*UR C RHO DEFINED RHO(1,1) = & -1 RHO(1,2) = & 0 C F2 DEFINED if (i8z.eq. 2) yd8z = & UIxx + 2*rk*URx - rk**2*UI - W*UI C RHO DEFINED RHO(2,1) = & 0 RHO(2,2) = & -1 call dtdpsx (yd8z,i8z,uu8z,un8z,rho,3,neqn) else endif endif return end function u8z(i8z,x,t0) implicit double precision (a-h,o-z) common/parm8z/ pi,p,rk ,a ,b u8z = dtdpxx() C############################################################################## C If you don't want to read the FINE PRINT, default the initial values. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + Now the initial values for the inverse power method may be defined +# C + using FORTRAN expressions. They may be functions of X. +# C + +# C + By default, the initial values are generated by a random number +# C + generator. This virtually eliminates any possibility of convergence +# C + to the wrong eigenvalue, due to an unlucky choice of initial values. +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C UR0 DEFINED C if (i8z.eq. 1) u8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] C UI0 DEFINED C if (i8z.eq. 2) u8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] return end subroutine gb8z(gd8z,ifac8z,i8z,j8z,x,t,uu8z) implicit double precision (a-h,o-z) parameter (neqnmx= 99) dimension uu8z(3,neqnmx) C un8z(1,I),un8z(2,I),... hold the (rarely used) values C of UI,UIx,... from the previous iteration or time step common /dtdp4x/ un8z(3,neqnmx) double precision none common/parm8z/ pi,p,rk ,a ,b none = dtdplx(2) UR = uu8z(1, 1) URx = uu8z(2, 1) UI = uu8z(1, 2) UIx = uu8z(2, 2) if (j8z.eq.0) gd8z = 0.0 C############################################################################## C Enter FORTRAN expressions to define the boundary condition functions, # C which may be functions of # C # C X,UR,URx,UI,UIx and (if applicable) T # C # C Recall that the boundary conditions have the form # C # C G1 = 0 # C G2 = 0 # C # C Enter NONE to indicate "no" boundary condition. # C # C +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++# C + If "no" boundary condition is specified, the corresponding PDE is +# C + enforced at a point just inside the boundary (exactly on the +# C + boundary, if EPS8Z is set to 0 in the main program). +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## if (ifac8z.eq. 1) then C############################################################################## C # C First define the boundary conditions at the point X = XGRID(1). # C############################################################################## if (j8z.eq.0) then C PERIODIC BOUNDARY CONDITIONS SET (SEE IPERDC) C G1 DEFINED C if (i8z.eq. 1) gd8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] C G2 DEFINED C if (i8z.eq. 2) gd8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] else endif endif if (ifac8z.eq. 2) then C############################################################################## C # C Now define the boundary conditions at the point X = XGRID(NXGRID). # C############################################################################## if (j8z.eq.0) then C PERIODIC BOUNDARY CONDITIONS SET (SEE IPERDC) C G1 DEFINED C if (i8z.eq. 1) gd8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] C G2 DEFINED C if (i8z.eq. 2) gd8z = C & [DEFAULT SELECTED, DEFINITION COMMENTED OUT] else endif endif return end subroutine pmod8z(x,t,uu8z,uprint,uxprint) implicit double precision (a-h,o-z) dimension uu8z(3,*),uprint(*),uxprint(*) common/dtdp14/sint(20),bint(20),slim8z(20),blim8z(20) common/parm8z/ pi,p,rk ,a ,b UR = uu8z(1, 1) URx = uu8z(2, 1) URxx= uu8z(3, 1) UI = uu8z(1, 2) UIx = uu8z(2, 2) UIxx= uu8z(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 UR,URx,UI,UIx at the output +# C + points. If different variables are to be saved (for later printing +# C + or plotting) the following functions can be used to re-define the +# C + output variables: +# C + define UPRINT(1) to replace UR +# C + UXPRINT(1) URx +# C + UPRINT(2) UI +# C + UXPRINT(2) UIx +# C + Each function may be a function of +# C + +# C + X,UR,URx,URxx,UI,UIx,UIxx and (if applicable) T +# C + +# C + Each may also be a function of the integral estimates SINT(1),..., +# C + BINT(1),... +# C + +# C + The default for each variable is no change, for example, UPRINT(1) +# C + defaults to UR. Enter FORTRAN expressions for each of the +# C + following functions (or default). +# C ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++# C############################################################################## C DEFINE UPRINT(*),UXPRINT(*) HERE: return end C dummy routines subroutine xy8z(i8z,iarc8z,s,x,y,s0,sf) implicit double precision (a-h,o-z) return end subroutine dis8z(x,y,ktri,triden,shape) implicit double precision (a-h,o-z) return end function fb8z(i8z,iarc8z,ktri,s,x,y,t) implicit double precision (a-h,o-z) fb8z = 0 return end function axis8z(i8z,x,y,z,ical8z) implicit double precision (a-h,o-z) axis8z = 0 return end subroutine tran8z(itrans,x,y,z) implicit double precision (a-h,o-z) return end subroutine postpr(tout,nsave,xout,nx,uout,neqn) implicit double precision (a-h,o-z) dimension xout(0:nx),tout(0:nsave) dimension uout(0:nx,2,neqn,0:nsave) common/parm8z/ pi,p,rk ,a ,b 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,IDER,IEQ,L) = U_IEQ, if IDER=1 C Ux_IEQ, if IDER=2 C (possibly as modified by UPRINT,..) C at the point XOUT(I) C at time/iteration TOUT(L). C ******* ADD POSTPROCESSING CODE HERE: C IN THE EXAMPLE BELOW, MATLAB PLOTFILES pde2d.m, C pde2d.rdm CREATED (REMOVE C! COMMENTS TO ACTIVATE) C! if (lun.eq.0) then C! lun = 46 C! open (lun,file='pde2d.m') C! open (lud,file='pde2d.rdm') C! 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!78754 format (i8) C! do 78755 i=0,nx C! write (lud,78760) xout(i) C!78755 continue C! do 78759 l=0,nsave0 C! write (lud,78760) tout(l) C! do 78758 ieq=1,neqn C! do 78757 ider=1,2 C! do 78756 i=0,nx C! write (lud,78760) uout(i,ider,ieq,l) C!78756 continue C!78757 continue C!78758 continue C!78759 continue C!78760 format (e16.8) C ******* WRITE pde2d.m C! call mtdp1dc(itype,lun) return end