LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAICN 510.84 T&63c vk>.6\-70. AUG. 31976 1 he person charging this material is re- sponsible for its return to the library from which it was withdrawn on or before the Latest Date stamped below. Theft, mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University. UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN XI a JAN. « 8 ^sx 5 Hfc t ; ?n '&# APR 8 I Wis , 3ft OCT EHOIO PHOTO REPRODLCTION OCT 2 2 SOT PHOTO REPRODUCTION NOV 1 6 fasdf REPRODUCE* 27KCI ftPTO0UCT(Oft sJeroouc 3CT 2 3 REtO L161 — O-1096 Digitized by the Internet Archive in 2012 with funding from University of Illinois Urbana-Champaign http://archive.org/details/useofilliacivtos67cich ENGINEERING LIBRARY UNIVERSITY OF ILLINOIS UR8ANA, ILLINOIS enter for Advanced Computation UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA. ILLINOIS 61801 CAC Document No. 67 THE USE OF ILLIAC IV TO SOLVE NUMERICAL FLUID DYNAMICS PROBLEMS By Paul T. Cichy February 27, 1973 CAC Document Wo. 67 THE USE OF ILLIAC IV TO SOLVE NUMERICAL FLUID DYNAMICS PROBLEMS A STATUS REPORT V Paul T. Cichv Center for Advanced Computation University of Illionis at Urbana-Champaign Urbana, Illinois 618OI February 27, 1973 This work was supported in part by the Advanced Research Projects Agency of the Department of Defense and was monitored by the U. S. Army Research Office-Durham under Contract No. DAHCOl*-72-C-0001, and bv the National Science Foundation (NSF GK 2813X). ABSTRACT The usefulness of ILLIAC IV in solving numerical fluid dynamics problems is studied by developing a parallel algorithm for the unsteady two dimensional flow of an incompressible fluid around a circular cylinder. Two cases are studied, the symmetric case, for which an operational GLYPNIR program is included, and the asymmetric case, for which a partially completed GLYPNIR program is attached. The programming of a parallel algorithm for ILLIAC IV is more difficult than programming a comparable serial algorithm but this additional difficulty is compensated for by the greatly reduced computation time required to solve the problem. TABLE OF CONTENTS Page 1. INTRODUCTION 1 2. THE ILLIAC IV SYSTEM 2 2.1 The ARPA Network 3 2.2 The Ames Site 3 2.3 The ILLIAC Array 6 2.k Programming Languages ° 3. THE CHOICE OF A PROBLEM 9 h. THE SYMMETRIC FLOW AROUND A CYLINDER 11 U.l The Problem 1] - ^.2 Development of the Algorithm 12 U . 3 Programming the Algorithm 24 k.k Testing the Algorithm 2U U.5 Problems with this Symmetric Formulation 2o k.6 Conclusions 27 5. THE ASYMMETRIC FLOW AROUND A CYLINDER 28 5.1 Description of the Problem 28 5.2 Development of the Algorithm 33 5«3 Programming the Algorithm 39 5.h Current Status of Program U8 5.5 Conclusions kQ 6. FUTURE PROGRAMS IN THIS PROJECT 1+9 REFERENCES 50 APPENDIX A. CICHY/GLYSOND 51 APPENDIX B. CICHY/GLYDIRECT 6h APPENDIX C. CICHY/MASTER 82 APPENDIX D. CICHY/DIRECT 105 APPENDIX E. CICHY/UNS 13^ 1. INTRODUCTION The numerical solution of a large class of fluid dynamic problems is currently impractical due to the enormous amounts of computation time needed to generate useful results on presently available serial processing computers. One approach to increasing computation speed using established hardware is to construct a computer which interprets a single instruction stream but causes arithmetic and logic operations to be performed on many sets of data simultaneously. This type of computer architecture has been given the descriptive terms "parallel processor" or "array processor". ILLIAC IV is a new parallel processing computer designed by the University of Illinois and Burroughs Corporation for the Advanced Research Projects Agency (ARPA). ILLIAC IV interprets a single instruction stream in its central control unit (CU) and then causes operations to take place in 6k separate processing elements (PE), each PE acting on its own set of data. The hardware was constructed by Burroughs at Paoli, Pennsylvania , and installed at the Ames Research Center, Moffet Field, California. The computer was accepted by ARPA in December 1972 and is expected to be operational by June 1973. The Ames facilities are accessable to a large group of users through the ARPA network, a national network connecting many Universities, Government Laboratories, and Government Contractors. There are several research areas where a very fast computer is currently needed in order to solve problems requiring the processing of large amounts of data. Problems in numerical fluid dynamics, studies in weather prediction, analysis of seismographic information to detect nuclear explosions, and the study of artificial intellegence all represent areas of active interest. The project discussed in this report is directed toward the solution of numerical fluid dynamics problems. The overall objectives are: 1) to determine what difficulties may be encountered in programming a parallel processing computer and to suggest ways to minimize them, 2) to develop a parallel algorithm which will solve a problem previously solved by serial techniques and then compare the two results for ease of programming, speed of obtaining the solution, and accuracy, and, finally, 3) to draw some conclusions about the general usefulness of ILLIAC IV in solving numerical fluid dynamic problems. Some of these objectives cannot be fully attained until ILLIAC IV is operational. This report presents the progress made toward meeting these objectives from June 1972 to February 1973. The initial section describes the ILLIAC IV architecture as it exists today. The following sections discuss the choice of a numerical fluid dynamic problem and the preparation. of two algorithms to solve that problem. These sections are followed by several conclusions on the use of ILLIAC IV based on the experience gained during this project . The project is jointly sponsored by the Department of Chemical Engineering (Dr. T. J. Hanratty) and the Center for Advanced Computation (Dr. A. Sameh). 2. THE ILLIAC IV SYSTEM The ILLIAC IV array processor does not exist as a totally separate entity but is intimately linked with peripheral equipment which manages its activity and performs the services it requests. The user approaches the ILLIAC IV system through the ARPA network and places a request for computing time with the TENEX managing system. TENEX in turn schedules and runs the job on the ILLIAC IV array processor and fulfills all the array's INPUT/OUTPUT (I/O) requests. A "brief description of each section of the ILLIAC IV system follows. A more complete description has been prepared by S. A. Denenberg [1]. 2.1 The ARPA Network The ARPA network is designed to make available to a large number of users a wide variety of computational resources. Some of the locations attached to the network are indicated in the sketch in Figure 1. Network users have available such large computers as IBM 360's, a Burrough's B67OO, and several PDP 10's and PDP 11 ' s . Each location is served by an Interface Message Processor (IMP) which controls the flow of information over the net. At present, the net is functional, but due to local conventions at each site, much care is needed in its use. A network protocal now under development should make its use much easier. 2.2 The Ames Site The Ames location is managed by a PDP 10 using a TENEX system as illustrated in Figure 2. Network users communicate with the TEN through a special control language. Documentation on this language is being prepared and should be available soon. The TENEX managing system controls a large disk data storage area, a UNICON laser storage device with a trillion bits of data storage and a B67OO for compiling programs for ILLIAC IV. It also manages the use of the ILLIAC IV array and handles all requests for I/O to and from the array. The TENEX managing system is being brought up slowly. At present it is capable of storing programs in its disk file Of o UJ CD o Q. HZ a: uj •— >- to < pqo: to < oa: ■^ o CNI UJ CD system and compiling ASK assembly language programs^but it is not able to run jobs on the ILLIAC array. Current projections are that it will be functional in April. Since all access to the ILLIAC array will be over the ARPA net and through the TENEX managing system, any use of the machine must wait until TENEX is functional. 2.3 The ILLIAC Array The architecture of ILLIAC IV is unique* and an understanding of this architecture is very important to those who would program ILLIAC efficiently. The computer is composed of two major sections as shown in Figure 3. The central control unit (CU) and the processing elements (PE), 6U units each having an arithmetic and logic unit (ALU) and a process element memory (PEM). The instruction stream is interpreted by the CU. The CU then sends messages to the processing elements to conduct the indicated data handling and arithmetic operations. Each PE executes the same command but operates on the specific data stored in its own processing element memory. Thus, when repetitive operations are to be performed on different but similar data sets, ILLIAC may be used to greatly increase the speed of computation by operating on 6k data sets at a time. It is projected that ILLIAC will be capable of executing 200 million instructions per second. One of the major limitations of the ILLIAC array is the small size of each PEM. The maximum capacity of each PEM is 20^8 words, that is approximately a 2K memory in each processing unit. For many problems this means that the storage cannot be core contained and that time consuming exchanges of data between the array and peripheral storage units must be carried out. ILLIAC operates with a 6U bit word allowing kQ bits of mantissa and 15 bits of exponent. This large word size will allow increased accuracy CO CQ O oc o o D UJ Q- v Z _l O oc I- z o o CO OS o J— O rH <\J m < oca: on oc -l <<<< ID UU OO £ <<<< o u < m Q O H < UJ UJ QHZ 2!rt_J CO vo UJ Q_ < CQ OS CO X CD OC zo PS 3 UJ oz a: o >-co < CO CC UJ a:o < o UJ Of >- •— • LJL, 11 The asymmetric case had been treated by D. C. Thoman and A. A. Szewczyk [5] and the results were available in their report. This problem was more difficult requiring memory swaps and special treatment of the boundaries and, therefore, was attacked only after the symmetric algorithm was completed. However, the asymmetric case, where vortex shedding occurs, is more realistic physically when considering flows with Reynolds numbers in excess of kO . The development of the algorithm for the symmetric case is presented in the next section. The asymmetric case is discussed in a later section, k. THE SYMMETRIC FLOW AROUND A CYLINDER U.l The Problem The unsteady two dimensional symmetric flow of an incompressible fluid around a circular cylinder was choosen as the first problem to be studied. This problem was investigated previously by J. S. Son [h] and his theoretical development and boundary conditions were adopted directly. The governing equations were developed from the two dimensional Navier-Stokes equations. These equations were differentiated in such a way that they could be combined to eliminate the pressure terms. The resulting equation along with the continuity equation formed a pair of coupled nonlinear equations describing the fluid flow. These equations were converted to a stream function and vorticity formulation and non- dimensionalized. The independent variables were then transformed to adjust the coordinate spacing so that the boundary layer region would be adequately covered with mesh points. The flow was assumed symmetric, thus the transverse velocity and the vorticity were assigned a value of zero along the centerline of 12 the flow. The fluid velocity at the outer edge of the flow field was taken to he parallel to the incoming flow and the velocity through the cylinder surface was taken as zero. The surface vorticity was estimated from the stream function distribution near the surface by expanding the vorticity function in terms of the stream function in the vicinity of the surface using the definition of vorticity. The pertinent equations developed in Son's thesis and the appropriate boundary conditions are given in Figure 5. k.2 Development of the Algorithm 1+. 2.1 Introduction The parallel algorithm developed to solve this problem was patterned after the serial algorithm used by Son. The general solution scheme is shown in Figure 6. The flow is assumed to start impulsively from rest and the initial stream function values are taken as those for the potential flow case. The finite difference form of the stream function equation is then used to adjust these values to be compatable with the boundary condi- tions. The initial vorticity in the flow field is assigned a value of zero. The surface vorticity at time zero is determined from the stream function distribution. When all the initial values have been determined, dimensionless time is advanced by one increment and a finite difference form of the vorticity transport equation is used to calculate vorticities at the new time level. The stream function values are then advanced by a finite difference form of the stream function equation, using the latest values of vorticity. Two different techniques for updating the stream function were used, an iterative method and a direct method. These will be discussed in more detail in a later section. 13 LU CD 11+ Q LU 1- >- H LU Z t: •— • ID o z O GJ t— 1 ►— • o ►— ■ LU a 1- U •—• 1- a: z ^?7 o •— • H cc LU n: z 1- o o z 1— 3 cc 3g > LU DC u. o 3 X O c2 > U. LU ^ o ^* 111 CJ a. CD LU __J LU »— • LU «^ LI- h- H UD _l t/> H or h- 0) LU «=c Z o to LU s: LU z: -I LU u_ LU ££ ■— • C£ CD < S LU 5 i- — * i— « •— • LU LU H H CD h- h- cc > < 0- z i- i— i «=c ►— • o -I Q s: •— * X U_ f— Z z o 0. o on LU Z3 i— i t— « to 3 o a. Z D_ t 1 O <_) 15 After the vorticity and stream function values have "been updated in the interior of the flow region, the surface vorticity is adjusted to he compatible with the current stream function distribution near the surface. If requested, the value of all variables at this time level can he printed out at this point. The computation can then proceed by incrementing time and repeating the steps described above or it may be terminated. The following sections treat in more detail each of the steps in the basic algorithm. k.2. 2 The initial solution The values of the stream function at time zero were obtained from the analytic equations for the potential flow of a fluid of infinite extent about a circular cylinder. Since the fluid is assumed to be set into motion impulsively, at time zero, the boundary layer will be confined to a vanishingly thin region next to the surface of the cylinder and will not affect the rest of the flow field. Since viscous effects are confined to the boundary layer, the motion of the majority of the flow field is inviscid and the potential flow solution describes the motion well. Due to the need to deal with a boundary which is not at infinity but at some large but finite distance from the cylinder, a slight adjustment must be made in the analytic infinite extent values of the stream function. This is accomplished with the use of the finite difference form of the stream function equation as described below. The actual surface vorticity at time zero should be infinite since, as the fluid impulsively assumes a velocity, a momentary velocity discontinuity exists at the surface. However, in order to obtain a tractable numerical solution, a large but finite vorticity is calculated from the stream function distribution. 16 The accuracy of the numerical solution at short times is in question due to the inability of the numerical scheme to account rigorously for the infinite surface vorticity and due to the difficulty in resolving the very thin boundary layer during these first few time steps. The effect that this initial error might cause in later calculations has not been determined but it is expected that it is confined to short times. One method to improve the accuracy of the short time solution would be to use an analytic representation of a growing boundary layer to deter- mine the values of the stream function and vorticity near the cylinder surface at a time when the boundary layer was still small with respect to the cylinder radius but when it was large enough to be resolved by the numerical scheme. The surface vorticity at this time would be finite. U.2. 3 The vorticity transport equation The Alternating Direction Implicit (ADl) method was used to obtain the values of vorticity at the next time level. This method has been shown to be stable for fairly large time increments. Two steps are required to reach the next time level. The finite difference equations for each step are shown in Figure J. During the first half time step the solution is taken to be implicit in the rows and explicit in the columns. In the next half time step the solution is implicit in the columns and explicit in the rows. This formulation results in the need to solve a tridiagonal matrix to obtain the implicit vorticity values in each row or column as the case may be. The Thomas Method represents a particularly simple method of solving these tridiagonal matricies. The solution is obtained by a series of calculations involving sequential terms in the coefficient matrix. The full power of ILLIAC IV is realized only when all 6k PE's are active 100% of the time. Although this is seldom accomplished, maximum IT C3 an. — • o -I q. a. CO S I— o i — i LJ O Q I— z QC U_ LU =3 O -J X O \- oo tu s < ~ az CD UJ 18 activity is a constant goal. Frequently the manner in which data is stored in PEM will determine how efficiently the computations can be carried out. For instance, solution of a tridiagonal matrix by the Thomas Method is inherently a serial process. When the solution is implicit in the columns the coefficient vectors for each column can be placed in a separate PE. Thus if the field width is choosen to be 62 (allowing two PE's for boundary conditions) then maximum use of the array processor can be obtained. The sequential processes take place as though 62 serial computers were calculating the result and thus 62 results are computed as fast as one. When the solution is implicit in the rows a problem arises. Since the process must follow a given sequence beginning with the first element in the row and continuing to the end of the row and then return, how can the first element of more than one row be accessed simultaneously? This problem is simply solved by shifting the coefficients in the second row one PE to the right, in the third row 2 PE's and so on until the entire field is skewed as shown in Figure 8. Then the first element in each of the 62 rows can be brought into separate PE's simultaneously. By appropriate indexing and routing of data among PE's, the solution for 62 rows is obtained simultaneously, thus using the array processor efficiently. h.2. k The stream function equation Two different techniques were used to update the values of the stream function, an iterative method and a direct method. Son's algorithm used only an iterative method. Timing studies showed that updating the stream function was the most time consuming step in the algorithm and thus 19 ^ Q 111 Ul CO < Q OO UJ CD m, _ — _ — V V V 1 * " ' T ' ' ' ' ft ■ ""•* ~ ■* J r f f ^s> CD »— « < a: \- co Q HI a: o »- CO < < Q V^ 20 prompted a search for a faster method. The direct method appears to cut the computation time 50 to 60 per cent. The iterative method used in Son's serial algorithm was the successive overrelaxation method (SOR). The computational molecule is shown in Figure 9-A and uses two updated values and two old values to determine the new center value. The serial algorithm thus starts in the lower left corner of the mesh and sweeps from left to right on each row, moving up the mesh. In order to use ILLIAC efficiently, a whole row must be processed at once so that the 62 PE's are all active. The first scheme tried used the computational molecule in Figure 9-B. Only the point from the previous row was an updated value. Unfortunately this method did not converge. A slightly more complicated but rapidly convergent method was devel- oped and will be called the modified successive overrelaxation method (MSOR) This method identifies each point in the mesh as either an odd point or an even point as determined by the sum of its indicies . The method requires two passes through the mesh. On the first, each odd point is updated using only old values as shown in the left computational molecule in Figure 9-C The second pass then updates the even points using all new information as shown in the molecule on the right. This staggering of points along with the proper index modification allows all the odd or even points in two rows to be updated simultaneously. The finite difference equation and a sketch of the computation mesh are given in figure 10-A. The MSOR method has been shown to have the same asympotic rate of convergence as the SOR method. A Center document by J. Ericksen [2] treats this in more detail. The MSOR technique makes efficient use of the ILLIAC IV array processor. 21 m*i ft U/' '-J A. POINT SOR /n srt m ti B. ROW SOR r V X (^ . X- OLP Ci MSOR /»»*/ m*i FIGURE 9 22 i '. \ FOURIER ANALYSIS FFT SOLVE ODE INVERSE TRANSFORM ^*"~ . MM^^BM^ ^ ^_ i- > \ J( i - O UI UJ z •-• o Q O CO C3 O to GO t— • o Q- K O _J CD GO CD UJ > •— ■ < a: UJ H o CO w 23 The stream function equation is Poisson's Equation and its iterative solution has received much attention. Recently Hockney [3] has proposed a technique for obtaining a solution by a noniterative method with a sig- nificant savings in computation time. Since the iterative solution of Poisson's Equation represented more than 80% of the total computation time for the problem, it was decided to try Hockney' s method. This direct method limits the number of mesh points in each direction to (2 +l), where n is any integer. This is an undesirable feature but in most cases it is quite tolerable given the increased computational speed. The program was initially developed in ALGOL and then converted to GLYPNIR . The GLYPNIR program was developed by J. Ericksen and is discussed in more detail in his report [2]. The method Fourier transforms the vorticity based variables along a whole row, resulting in a series of Fourier coefficients which are independent from column to column (PE to PE) but are implicit up the columns. The computational mesh is shown in Figure 10-B. In effect the Fourier analysis has reduced the problem from one of solving a partial differential equation to one of solving a series of coupled ordinary differential equations, This later solution is accomplished by inverting a simple tridiagonal matrix. The resulting inverse Fourier components are then converted back by Fourier synthesis to yield the stream function distribution. The method reportedly uses ILLIAC efficiently and thus should represent a savings of time over the iterative method. Test programs for each method of updating the stream function have been prepared in GLYPNIR . The iterative program is titled CICHY/GLYSOND and the direct program is called CICHY/GLYDIRECT. 2U h. 3 Programming the Algorithm The strategy used in programming this algorithm was to use subroutines to represent separate logical steps in the algorithm. The subroutine MAIN controls the program, calling other routines as needed. The subroutine names are descriptive of the function they perform. The program is designed to handle up to 62 internal mesh points across the field and up to 25^ internal mesh points along the field. Boundary conditions are inserted outside these points to complete the definition of the field. The program parameters are specified in the subroutine SYSTEM PARAMETERS. They are the Reynolds Number (NREY), the dimensionless time increment (DT), the dimensionless time limit for the problem (TLMT), the error limit for the convergence of the iterative solution of the stream function equation (LIMIT), the total (includes boundary) number of mesh points in the circumferential direction (NETA), the total number of mesh points in the radial direction (NTSl), the maximum value of the circumferential coordinate (ETAMAX), the maximum value of the radial coordinate (TSIMAX), and the increment in time steps at which printouts are to be taken (PINC). Given these parameters > the program generates all the other information and control parameters it needs to solve the problem. The two GLYPNIR programs developed are listed in Appendicies A and B. These programs contain many comments indicating how the actual operations of the algorithm are carried out in the computer. ^.^ Testing the Algorithm GLYSOND is coded in GLYPNIR and has been successfully compiled on the GLYPNIR Compiler. The ASK produced from the compilation has been 25 successfully compiled into machine language. This machine language has not "been tested on the ILLIAC IV. However, a simulator for the ILLIAC IV hardware has been developed and runs on a Burroughs B67OO. With this simulator, machine code can be tested just as if it were running on ILLIAC IV. The major problem with testing code on the simulator is that, since the simulator runs in a serial manner, it runs much slower than ILLIAC , in fact the ratio is 200,000 to one. It was not possible to test the machine code from GLYS0ND for a reasonably sized mesh. Therefore, a 9x17 mesh was chosen and the number of iterations in the iterative MS0R was limited to six. Even so, to accomplish three time steps on the simulator required 12 minutes of CPU time. In order to check the results of this abbreviated mesh problem it "Was necessary to develop a serial solution to the problem. Son's FORTRAN program was converted to an ALGOL program titled CICHY/MASTER, Appendix C. MASTER was used on a 31x51 mesh at Reynolds Number of U0 and the results checked well with those of Son. The same program was then used to generate results for 3 time steps using a 9x17 mesh and with iterations in the S0R limited to 6. These results compared well with the results from GLYS0ND. This close comparison was taken as proof that the GLYPNIR program worked well and was ready to be tested on ILLIAC IV. The same test procedure was applied to the program using the direct solution of the stream function equation. The GLYPNIR program is titled CICHY/GLYDIRECT, Appendix B, and the ALGOL program is CICHY/DIRECT, Appendix D. In order to have numerical results to check the solution obtained from GLYDIRECT on ILLIAC IV, DIRECT was used to solve the problem on 26 a 33x65 node mesh. These numerical results should compare closely with the numerical results obtained from GLYDIRECT run on ILLIAC IV using a similar sized mesh if the program is operating correctly. These numerical results are currently stored with the numerical results of Son at the Department of Chemical Engineering. The next step in testing GLYSOND and GLYDIRECT is to run them on ILLIAC IV. This should be possible at least by June 1973- GLYSOND as shown in Appendix A is set up to solve a 9x17 mesh and limit MSOR to six iterations. A copy of this program is now resident in the TENEX file system at Ames and John Gaffney and Marvin Graham will attempt to run this as soon as the ILLIAC is ready. These results should compare well with those obtained from the simulator (the printouts are stored at the Department of Chemical Engineering). The next step should be to remove cards 153 and 597 from GLYSOND which will allow MSOR to iterate to con- vergence. Then run GLYSOND with NREY=l+0, NETA=31 and NTSI=51, and compare these results with those obtained by Son. CPU time should be requested so that comparisons of running time and accuracy can be made. After GLYSOND has been tested, GLYDIRECT should receive similar testing. Because of changes in the GLYPNIR compiler one or two variable names in GLYDIRECT may now be reserved words. The compiler will point these out and they should be changed. At present the final 1/0 specification are still not complete. It is expected that they should be finished soon and that as soon as they are finished a form of GLYPNIR 1/0 for ILLIAC IV will be implemented. h . 5 Problems with this Symmetric Formulation There are several features of the symmetric formulation of the problem which indicate that a more general approach is desirable. 27 The case of symmetric vorticies is realized in practice only up to a Reynolds number of about kO. Above this Reynolds number the flow breaks down into a pattern of vortex shedding. The symmetric formulation cannot represent this flow pattern and is thus limited in its range of application. The use of circular cylindrical coordinates limits the extent to which the downstream boundary can be moved away from the cylinder. The parallel flow assumption made at this boundary is fairly restrictive and it is important that it have a negligible effect on the size and shape of the trailing vorticies and on the pressure distribution around the cylinder. Due to the diverging radial coordinate lines, an unreasonable number of mesh points must be added to retain adequate downstream resolution as the boundary is moved away from the cylinder. The mesh spacing in this formulation follows a strict geometric pattern, modified slightly to better account for the boundary layer. In fact, there are regions other than the boundary layer where variable gra- dients may be large but which are not adequately covered with mesh points, such as the shear layer at the near edges of the vorticies. It would be desirable to be able to tailor the mesh size in any region to the size of the gradients expected in that region. This control over mesh size would require a significantly more complicated program. h.6 Conclusions The ability of a parallel processor to solve one particular type of numerical fluid dynamics problem has been established with the completion and testing of the GLYSOND and GLYDIRECT parallel programs. With the proper choice of the number of mesh points, these algorithms should operate effi- ciently and result in very rapid solutions. Actual testing on ILLIAC IV will be required to determine just how much faster the parallel program is than the serial program. 28 The placement of data in the PEM was shown to be very important in maintaining a high efficiency during certain portions of the program. The routing overhead required is generally negligible, except possibly in the case of very short rows. The programming of a parallel machine is more difficult than a serial machine, but given a good understanding of the machine architecture and a working knowledge of both programming languages, the programming of an algorithm such as this one is not unreasonably difficult. Planning of data storage areas, enabling mode control, and data movement patterns represent some of the more difficult tasks. Some of the most difficult problems encountered in the development of the symmetric algorithm related to the uncertainties in the supporting software and the unreliability of local hardware. Hopefully both of these situations have been improved to the point where they are no longer serious problems . 5. THE ASYMMETRIC FLOW AROUND A CYLINDER 3.1 Description of the Problem The asymmetric flow around a cylinder is of particular interest since flows of this type are frequently observed in actual fluid flow situations. The famous von Karmen vortex street is an example of asymmetric vortex shed- ding behind a cylinder. One objective of this second study, besides relaxing the conditions of symmetry, was to find a coordinate system which would allow the downstream boundary condition to be moved away from the cylinder while not increasing unreasonably the number of mesh points needed to maintain resolution. This constraint suggested that some kind of a rectangular coordinate system 29 would "be appropriate downstream. However, near the surface of the cylinder, where resolution is needed for the boundary layer, a circular cylindrical coordinate system seemed best. The search for a coordinate system which would be cylinder like near the cylinder surface and rectangular downstream led to the trial of a set of coordinates based on the potential flow solution. These conformed well to the shape of the cylinder and became rectangular at large distances from the cylinder surface. These seemingly perfect coordinates had, however, one fatal flaw, a singularity at the leading and trailing stagnation points. After considering several methods to circumvent the singularity these coordinates were abandoned in favor of the combination of circular and rectangular coordinates proposed by D. C. Thoman and A. A. Szewczyk [5]. The theoretical development used for the asymmetric problem has been taken directly from the report mentioned above. The equations to be solved are basically the same as those developed in Son's work but the boundary conditions are slightly different since the flow is asym- metric. The boundary conditions and the organization of the flow field are shown in Figure 11. The five regions shown in Figure 11 represent areas of differing fluid activity. The mesh distribution in each region has been carefully chosen to adequately resolve the expected variable gradients. The mesh distribution over the important parts of Regions 2 and 3 is illustrated in Figures 12 and 13. Although these drawings are not exact, they show the increased concentration of mesh points just behind the cylinder, especially in the region where a shear layer might be expected. Fewer points are placed in Region 2, upstream of the cylinder, since the flow is expected to be quite regular there. 30 CO o M Eh H P O O K O PQ P CO o H O « J* o Fn > M 31 5.0 4.0 3.0 - 2.0 1.0- 0.0 + 1 - 1 1 1 1 1 J 1 1 — 1 — I 1 — I — J-J — I 1 1 ! A 1 1 , -5.0 -4.0 -30 -2.0 -1.0 0.0 PART OF REGION 2 FIGURE 12 32 5.0 4.0 3.0 2.0 1.0 0.0 , -1.0 • nz S9BHB is:: :':: ITT III -Ij -Ij E::_: z.l t~: : III t" . ~.t EZ- " -It t;:z: ::± tz: _ : OJ0 1.0 2D 3.0 4.0 PART OF REGION 3 FIGURE 13 33 The relative sizes of the lateral mesh spacing used in Regions 3, h and 5 are shown in Figure lk . The length of Regions k and 5 can be changed easily to assess the effect of the downstream boundary condition on the solution. The mesh spacings used in the cylindrical coordinates of Region 1 are illustrated in Figure 15. The points at R=2 correspond to points in the rectangular meshes of Regions 2 and 3. 5.2 Development of the Algorithm The general algorithm for the solution of the asymmetric problem follows the same pattern as discussed for the symmetric problem. A list of the steps in the algorithm are given in Figure l6. Two sets of finite difference equations were required, one in rectan- gular coordinates and the other in cylindrical coordinates. The finite difference forms of the vorticity transport equation and the stream function equation for rectangular coordinates are shown in Figure IT. The equations are explicit and allow for a variable mesh size through the inclusions of the cell width b. and cell height a., i J The presence of the circular coordinate system in the interior of the rectangular system excludes the use of implicit techniques. It is, therefore, necessary to devise a finite difference form of the vorticity transport equation which is explicit. This form requires that two com- ponents of the fluid velocity, u and v, be computed from the stream function distribution before the vorticity values at the next time level can be obtained. The implicit methods used to solve for the vorticity in the symmetric problem were stable for almost any reasonable time step. The stable time step for the explicit method is generally smaller, thus increasing the computation time required to obtain a given solution. 3h 5.0 Y 3 Y 4 Y 5 4.0 3.0 ? 1.0 nn REGION 3 REGION U FIGURE Ik REGION 5 35 20 1.95 1.90 1.85. 1.80 1.75 1.70 1.65 1.60 1.55 1.50 1.45 1.00 135 1.30 1.25 1.20 1.15 110 1.05 1.00 0.0 2.5 .7 5 1.0 REGION 1 FIGURE 15 36 CO Ul CD Q ea UJ z Ul -J _l >- or Ul 1- o z u o z o a o_ •— « o \- •— « Ul UJ s CO *— « < H a: h- UJ H or ►— • co z 3 u CO o 13 3 _J z UJ > a Ul y^ u. < 13 •-* UJ t: 3C > u. V- UI oc —* H 1— s: UJ t— o i- ¥— « < s >- 5" o < z Q^ HI •— i 1- —t UJ -J ar X X o H ^ 1- O or Ul =3 ^ UJ 1 co i- ■— • f- > CO z Ul -J < UJ o > Ul 5 3 o 1- S ^^ 1- or < a. Q. z a. o •— • o ^ a s: s: »— • s: »-^ z z Ul Q- o o or o h- *— • »—* z r> o o fl- u 5 t 1 Q_ m± * CD C_> 37 o U_ UJ o z Ul uj IX. <3 Hi o z < cc I • -» V £ + I <3 -tf * :i o >- ^ i i ^ ^ VC" o > o Hi O z HI OH HI Q z % ~7 4. ^? 1 -JA * + -$ -* \ +• 4. » V?H ^ a HI z o =3 h. H i 38 Thoman and Szewczyk discuss the stability requirements for the explicit method in their report. The difference approximation to a partial differential operator on a nonuniform mesh has a larger truncation error than that on a uniform mesh of similar cell size. The additional error arises due to the fact that the coefficient of the first neglected term in the Taylors series approximation to the derivative contains the difference of the mesh sizes on each side of the point of interest. When the mesh size is uniform, this difference is zero and the effective accuracy of the approximation is increased "by one order. However, in the case of a nonuniform mesh, this term must be retained and it contributes to the error in the approximation. The importance of the error conditions mentioned above can be minimized when taken over the whole mesh by keeping the cell sizes constant wherever possible. When moving from a region of one cell size to that of another, the amount of cell to cell variation should be limited to an amount such that the coefficient of the first neglected term is always nearly zero, in addition, if the regions of changing cell size can be placed in areas where the value of the first neglected derivative itself is expected to be small, the error encountered will be even smaller. Another technique to reduce error has been proposed but was not used in this study. After applying the explicit form to compute the values at the next time step, use the finite difference forms of the neglected deriva- tives to evaluate the error caused by truncating the Taylors series approxi- mation. Subtract this error from the computed values and repeat the process until the calculated truncation error is reduced to zero. This technique would improve the accuracy of the symmetric and asymmetric solutions and consideration should be given to applying it to these problems. 39 5. 3 Programming the Algorithm The program for the asymmetric flow problem is coded in the latest version of GLYPNIR . Appendix E contains the sections of code which have been completed at this time. Further work on the program will be needed before it can be run on ILLIAC IV. The following sections of the report discuss some of the problems encountered in programming this large code for ILLIAC , and describes some possible solutions. No attempt to optimize the code has been made at this point in its development. 5.3.1 The buffer system The symmetric problem was small enough to be core contained. The asymmetric case, however, was too large to fit in the PEM and thus a buffer system was developed. The total number of rows of PEM storage required to hold the entire stream function mesh was first determined. The total was then divided into four sections each containing nearly the same number of rows. This division was made strictly on the basis of the number of rows, without regard to the boundaries of the various flow regions. The section sizes were then doubled to account for the vorticity storage needed. A separate disk file was set up to hold each of the k sections. Two areas were reserved in the PEM, each area large enough to hold one section. The Buffer System is illustrated in Figure 18. During the course of the computation, a section is read into or written out of one area while calculations proceed in the other. When calculations and transfers are complete, calculation moves to the second area and transfers begin on the area just calculated. Because of ILLIAC 's great speed, it is expected that the transfer times will exceed the calculation time in most problems, a condition known as being I/O bound. ko CO H w s H 1+1 A working row (or rows) was designated outside of the transferable areas in PEM to assist in the transportation of information from one section to the next. When calculations are complete in a section, the last row in that section is moved to the work area. The first operation on the next section is to insert the new information from the work area into the first row. The calculations can then proceed as usual. The sectioning of the program and the use of two buffers suggested that a separate piece of code could be used to carry out the calculations in each section. This idea has been implemented through the use of a separate subroutine for each section. The technique provides a logical organization to the program. It is particularly useful in this problem as it allows a fairly complicated program to be broken down into a series less complicated units. 3-3.2 Organization of PEM The arrangement of the mesh points in PEM is a very important operation in setting up the program. It is necessary to lay out each buffer as a separate unit, indexed all the way through with respect to its first row. Figures 19 through 22 present the allocation of PEM storage for each buffer load. The numbers along the right edge are the buffer indicies. On their right are the coordinate value of the mesh line stored in that PEM row. The numbers on the left represent the regional indicies, indexed relative to the first row in the region. All of these indicies are used in the program. The number of data points in a mesh row exceeds 6k in several regions. This makes necessary the sectioning of the row into pieces, 1+2 FIGURE 19 k3 FIGURE 20 kk /37\ 129 /SI •»"— ^mm^t i . i 1 i • i . ...^ j f i 1 i I I ;■ I I! 23* ■> 2**7 ZHZ 2VJ } 23.* i/M 3l i i ■ ■ I i i i i__i. 4. FIGURE 21 1*5 * FIGURE 22 he each of which will fit into a PE row. This sectioning results in complicated indexing problems. The sketches shown in Figures 19 to 22 are very useful in keeping the indexing straight. It is also necessary to leave at least one working position each side of the split in the row and to exchange data between working spaces as the computation proceeds from one PE row to another within a mesh row. This does cut into the efficiency of the calculations but it is difficult to assess how much. The boolean mode variables control the enabling and disabling of PE's in the array. Each row has its own mode pattern. These patterns are carefully constructed so as to define exactly the active mesh. Separate mode patterns must be used when boundary conditions are evaluated. The data allocation sketches are very helpful in setting up the proper mode control patterns. The nontransferable portion of PEM contains a number of variables which are associated with specific mesh points. The principle variables are the coordinate values of each mesh point, the cell width and height of each mesh point and the modes for each row. These are all assessed by regional indicies and are given a variable name which indicates to which region they belong. The coordinates and cell sizes must be provided by the user before the program is compiled. A separate procedure could be developed to obtain the cell sizes from the coordinates but that is not implemented in this program. 5. 3. 3 Boundary conditions The boundary conditions which involve the endpoints of PEM rows require special attention. Each point must be treated separately, as in a serial machine. Since the endpoints are a small fraction of the total hi mesh points, this inefficiency may not be too serious. The "boundaries which involve whole PE rows are handled with little loss of efficiency. The downstream "boundary conditions used in this problem were taken directly from the report of Thoman and Szewczyk. They are that the axial velocity gradient and the axial vorticity gradient are both zero. These boundary conditions were selected from those experimentally tested as the ones which produced the most desirable solution. These conditions should be used in the initial version of this program but other more physically realistic conditions should be examined in the future. A boundary layer type of outflow condition with — ^ = and — — = might 9: d y. produce satisfactory results. The fact that the distance between their cylinder and the downstream boundary was only 8.5 diameters may be one reason why Thoman and Szewczyk had difficulty with other boundary conditions. The asymmetric program developed here extends to 50 diameters. 5.3.^ Changing regions At the place where the mesh distribution across the rows changes, a special treatment of the data must take place. If more mesh points are added, an interpolation is required. This is the case in going from Region 2 to Region 3. It can be seen that the new points needed have been placed exactly halfway between the existing points so that a simple averaging can be used to obtain the new values. The process is as follows. The last row in Region 2 is moved to a work array. A second work array is filled consecutively with values from the first array but every other position is left blank. From the values on each side of a blank an average is obtained and inserted in the blank space. This new work area is then moved to the first row in Region 3 and the calcu- lations proceed as usual. hQ This method works even when simple averaging is not used. The appropriate number of blanks must be left in the new work area and inter- polation on the bracketing values must then be carried out using the coordinate values of the points involved. 5.1* Current Status of Program The overall computational scheme to be used in the solution of the asymmetric problem has been developed and can be seen in the structure of the procedures used in CICHY/UNS in Appendix E. Much of the detail discussed in the preceding sections has yet to be programmed. Keeping track of indicies and control modes represents a most difficult task but frequent reference to the allocation sketches has proved very helpful. 5 . 5 Conclusions The asymmetric problem is an excellent test case for evaluating the use of ILLIAC IV in solving fluid dynamic problems. The features which make it most attractive are its large size, requiring the use of buffers and overlapping 1/0, and its nonuniform mesh size, with the accompanying indexing problems. These features are held in common with many of the types of fluid dynamic problems invisioned for solution on ILLIAC IV. Thus if this problem can be solved efficiently and rapidly, there is great hope for the increased use of ILLIAC IV in numerical fluid dynamics . The programming of a large problem such as this can be speeded by the use of storage allocation charts. These are particularly helpful in determining the correct indicies and enabling mode patterns. The use of a buffer system to handle problems with large data sets appears to work well. Only when actual 1/0 times and processor times are available can the actual efficiency of this method be determined. h9 6. FUTURE PROGRAMS IN THIS PROJECT There are several interesting studies which can be "built on the present work. First, of course, the asymmetric program should be completed. This, along with the symmetric program, should then be run on ILLIAC IV and careful account made of the time actually needed to solve the problem. The resulting solutions should then be compared to the serial solutions previously obtained and some conclusions drawn as to the benefit in speed and accuracy realized from using ILLIAC IV. The effect of the actual mesh point density and distribution could be carried out easily once the asymmetric program is completed. It would be interesting to determine how to construct a mesh with just enough points to resolve the fluid motion accurately. The effect of several downstream boundary conditions on the numerical solution could also be evaluated using the asymmetric program. It would be particularly interesting to evaluate boundary layer type conditions, since they imply that any fluid motions beyond the downstream boundary will not feed back to affect the solution near the cylinder. An entirely different numerical approach to the problem of flow around a cylinder might be tried. Recently finite element methods have emerged as possible ways to solve fluid flow problems. Adapting these techniques for use on ILLIAC IV could be a useful and rewarding endeavor. 50 REFERENCES [l] Denenberg, S. A., "An Introductory Description of the ILLIAC IV System," CAC Document No. 10, Center for Advanced Com- putation, University of Illinois at Urbana-Champaign. [2] Ericksen, J., "iterative and Direct Methods for Solving Poisson's Equation and their Adaptability to ILLIAC IV", CAC Document No. 60, (1972), Center for Advanced Computation, University of Illinois at Urbana-Champaign. [3] Hockney, R. W., "The Potential Calculation and Some Applications", Methods in Computational Physics, Volume 9« [k] Son, J. S. and Hanratty, T. J., "Numerical Solution for the Flow Around a Cylinder at Reynolds Numbers of U0, 200 and 500", J. of Fluid Mechanics 35 (2), 369 (1969). [5] Thoman, D. C. and Szevczyk, A. A., "Numerical Solutions of Time Dependent Two Dimensional Flow of a Viscons, Incompressible Fluid Over Stationary and Rotating Cylinders", Technical Report #66-lU (1966), Heat Transfer and Fluid Mechanics Laboratory, Department of Mechanical Engineering, University of Notre Dame. 51 APPENDIX A $2IP 00000100 JLIBRARY 00000200 * CICHY/GLYSOND 00000300 % 5/23/72 00000400 BEGIN 00000500 CREAL Pit DETA. DTSI. DT. OX. IOX. RLXP. LIMIT* NREY. KE. 00000600 TLMT. ORE* ETAMAX. TSIMAX, TERM4I .TERMB. TERMB2. 00000700 TERMD.TERMD2.0TI2.AT.BT.FT.CT.ET.GT.TA.DXS; 00000800 CINT NETA. NTSI. LL1. LL2.LL3.LL4.LL5.LL6. NITER. PINC.LM . 00000900 EVEN I 00001000 PINT CEVEN. CODD; 00001100 PREAL ETA. M. S ; 00001200 PREAL VECTOR PSII 201. ZETAl 201 I 00001300 BOOLEAN BCS1. BC11. BCS. SAVEMODE. BC» AD.BCT. BEV.BOD.BEVEN. 00001400 BODD* 00001500 FILE LINE (10 ROWS ) * 00001600 FILE I4DISK1 =»A / DISKFILE'M 90 ROWS FULL ) I 00001700 LABEL FINI I 00001800 % 00001900 % 00002000 % 00002100 % ««««««««««»«»««««»«««««««««»«««««»«««««>«>«>««««»««>««««»«««>»««««»«»««««00002?00 % 00002300 PE REAL SUBROUTINE 00002400 EXP AS RGA (PE REAL ARG AS RGA) I 00002500 % EXPONENTIAL 00002600 BEGIN 00002700 $* CALL EXP64 I 00002A00 END » 00002900 % •«««»««*»«««««*«««««««**«««*«««««**»»««»««««««»««»»««««*«*««««««««««00003000 % 00003100 PE REAL SUBROUTINE 0000320C SIN AS RGA (PE REAL ARG AS RGA)) 00003300 % SINE 00003400 BEGIN 00003500 $« CALL SIN64 I 00003600 END t 00003700 % «»»«»»««»««««««««««*»««««««»«»««*»«*««*««««««*«»«»*««««««*«««*»»»*««00 003800 % 0000390C PE REAL SUBROUTINE 00004000 LN AS RGA ( PE REAL ARG AS RGA ) I 00004100 % NATURAL LOGARITHM 00004200 BEGIN 00004300 $* CALL LN64 I 00004400 END ; 00004500 % *«<»««k»«»«»»»«> < h»»<>* »»•*♦» »..u •», •••».. »<>.••««•••••••••«•«••••» 004600 % 00004700 PE REAL SUBROUTINE 00004800 ARCTAN AS RGA (PE REAL ARG AS RGA) I 00004900 % ARCTANGENT 00005000 BEGIN 00005100 $» CALL ARCTAN64 ; 00005200 END ? 00005300 % ««««»««»««»»««*««««« «««»««u»««*«»«««««*«»« »««••»«««««•»«»••«*••«*»«« 00 005400 % 00005500 PE REAL SUBROUTINE 00005600 PRINT1 (PCDOINT VAR. CINT PRNT) ; 00005700 % PRINTS VALUES OF VARIABLE VAR FOR 00005800 % THE ENTIRE NET 00005900 BEGIN 00006000 52 CINT KL ; % PARAMETERS REQUIRED ZETAIKL1.LL6 SAVEMODE := MODE ; MODE := BCS ; LOOP KL := 0, 1, LL6 DO BEGIN IF PRNT EQL 1 THEN SIMWRITE (LINE* "VORTICITY* ROW NUMBER "♦ KL ) ; IF PRNT EQL 2 THEN SIMWRITE (LINE* ••STREAMFUNCTION* ROW NUMBER ••♦ KL ) f IF PRNT EQL 3 THEN SIMWRI TE (LINE » ••INITIAL STREAMFUNCTION* ROW NUMBER "tKL) J SIMWRITE (LINE* VARIKL) ) ; END ; MODE := SAVEMODE ; END ; % PE REAL SUBROUTINE SYSTEMPARAMETERS ♦ BEGIN MODE NREY DT : = TLMT LIMIT NETA NTSI := TRUE := 40.0 0.04 := 0.08S := 0.00001 ; := 9 ; % NUMBER OF THETA MESH POINTS := 17 ; % NUMBER OF RADIAL MESH POINTS a := l.o ; x := 1.5 ; := 1 ; ETAMA TSIMA PINC END ; COMPUTESYS BEGIN MODE PI : BCS1 BC11 AD : DETA DTSI LL1 PE REAL SUBROUTINE TEMCONSTANTS ; := TRUE ; = W » ARCTAN(l) ; := BOOLEAN(8000000000000000(16) ) ; := BOOLEAN(0C000000000000000(16) ) 5 = BOOLEAN(0000000000000001 (16) ) ; := ETAMAX /(NETA - 1) ; := TSIMAX /(NTSI - 1) ; := NTSI - 2; % LOOP LIMIT ONE* THE NUMBER OF * INTERNAL * % MESH POINTS IN THE RADIAL (TSI) D1RECTI LL2 := NTSI - 3; % LOOP LIMIT TWO* INDICATES TO THE PROGRAM % WHEN THE LAST SECTION HAS BEEN PROCESSED* % OR* WHEN THE LAST ROW HAS BEEN PROCESSED % IN THE THOMAS METHOD. IF LL1 LEQ 62 THEN LL3 := LL1 ELSE BEGIN IF LL1 LEQ 124 THEN LL3 := LL1/2 ELSE BEGIN IF LL1 LEQ 186 THEN LL3 := LL1/3 ELSE BEGIN 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 *»000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 **ooo 000 000 000 000 000 000 000 000 000 000 000 000 000 ON000 000 000 000 000 000 000 000 000 000 000 000 06100 06200 06300 06400 06500 06600 06700 06800 06900 07000 07100 07200 07300 07400 07500 07600 07700 07800 07900 08000 08100 08200 08300 08400 08500 08600 08700 08800 08900 09000 09100 09200 09300 09400 09500 09600 09700 09800 09900 10000 10100 10200 10300 10400 10500 10600 10700 10800 10900 11000 11100 11200 11300 11400 11500 11600 11700 11800 11900 12000 53 IF LL1 LEQ 248 THEN LL3 := LL1/4 ELSE BEGIN GO TO FINI END END END END J % LOOP LIMIT THREE* THE NUMBER OF INTERNAL % RADIAL MESH POINTS IN A SECTION* GIVEN % HERE AT ITS MAXIMUM VALUE OF 62. THIS % NUMBER MUST BE DETERMINED BY THE USER T % FIT THE NET LENGTH SO THAT LL1/LL3 IS % A WHOLE NUMBER. LL4 := NETA - 2 LL5 := NETA - % LOOP LIMIT FOUR* INDICATES THE NUMBER OF % INTERNAL MESH POINTS IN THE CIRCUMFEREN % (ETA) DIRECTION. THE MAXIMUM NUMBER IS 3 I % LOOP LIMIT FIVE* INDICATES WHEN THE LAST % COLUMN HAS BEEN PROCESSED IN THE % THOMAS METHOD LL6 := NTSI - 1 KE := PI » DTSI DX := EXP(KE) I IDX := 1.0/DX ; % % CALCULATE THE RELAXATION PARAMETER % BEGIN CREAL AR* BR* CR* OR I / ( NETA ♦ NETA ) / ( NTSI • NTSI ) ♦ BR ) • 0.5 I * LN(CR) ; RLXP := 1.84 = 1.0 = 1.0 = (AR = 0.5 := EXP(DR) ; RLXP := 2.0 / (1 ; end ; AR BR CR DR DRE ♦ PI • DRE > I % % CALCULATE THE BOOLEAN MASKS NEEDED FOR THE ARRAY WIDTH USED. 3EGIN CINT J ; 3C := FALSE ; LOOP j:xl.l»LL4 DO BEGIN BC := BC OR BCS1 BC := REVR(1*BC) END ; = REVR(1*BC) ; :■ BCS OR BC11 ; BCS BCS end ; % THE BOOLEAN VARIABLE BC MUST REFLECT THE % WIDTH OF THE NET. A "1" SHOULD BE PLAC % IN EACH BIT POSITION WHICH CORRESPONDS % TO A PE CONTAINING ONE OF THE INTERNAL % POINTS. ALL OTHER BITS MUST BE ZERO. 3EGIN CINT j; BCT := FALSE I LOOP J := 1.1*LL3 DO BEGIN BCT := BCT OR BCS1 i 000 000 000 000 000 000 000 000 000 000 000 TIAL000 62 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 ED 000 000 MESH000 uoo 000 000 000 000 000 000 12100 12200 12300 12400 12500 12600 12700 12800 12900 13000 13100 13200 13300 13400 13500 13600 13700 13800 13900 14000 14100 14200 14300 14400 14500 14600 14700 14800 14900 15000 15100 15200 15300 15400 15500 15600 15700 15800 15900 16000 16100 16200 16300 16400 16500 16600 16700 16800 16900 17000 17100 17200 17300 17400 17500 17600 17700 17800 17900 18000 5h BCT := REVR(ltBCT) ; end; % % end; TERM4I TERMB TERMB2 TERMD TERMD2 DTI2 = 0.25 /(DTSI » DETA) ; = 2.0 /(NREY * DETA * DETA) ; = 2.0 * TERMB ; = 2.0 /(NREY • DTSI • DTSI) ; = 2.0 * TERMD ; = 2.0/DT ; BEVEN := B00LEAN(2AAAAAAAAAAAAAAA(16) ) ; BODD := BOOLEAN (5555555555555554 (16) ) ; IF ( ( LL1 DIV 2) * 2 - LL1 ) = THEN BEGIN EVEN := 1 ; LM := LL1 ; END ELSE BEGIN EVEN is o ; LM := LL2 ; END ; BOD := BODD AND BC ; BEV := BEVEN AND BC ; CEVEN := ; CODD := ; dxs := dx » dx ; MODE := BEV : ceven := l ; MODE := TRUE ; MODE := BOD ; CODD := 1 ; AT := RLXP * DTSI • DTSI • BT := 2.0 » (DTSI * DTSI + CT := AT/BT ; FT := RLXP • DETA • DETA ; ET := FT/BT ; GT := CT * DETA » DETA ; TA := 1.0 - RLXP ; DETA » DETA) ; ETA MODE MODE := TRUE ; MODE := (PEN LEQ (NETA - 1) / 2 ) ; ETA := PEN ; MODE := TRUE ; MODE := (PEN GTR (NETA - 1) / 2 ) ; = (NETA - 1 ) - PEN ; :- BC ; M := ETA • PI * DETA ; S := SIN(M) ; MODE := TRUE ? END ; % PE REAL SUBROUTINE COMPUTEADICOEFFICIENTSANDSKEW (CREAL Et CNPOIN CNPOINT COFA» PCPOINT KFA, PCPOINT KFC* PCPOINT KFD ) BEGIN % SKEW COEFFICIENTS CINT K»J»L ; PREAL TERMA* TERMC ; 00018100 00018200 00018300 00018400 00018500 00018600 00018700 00018800 00018900 00019000 00019100 00019200 00019300 00019400 00019500 00019600 00019700 00019800 00019900 00020000 00020100 00020200 00020300 00020400 00020500 00020600 00020700 00020800 00020900 00021000 00021100 00021200 00021300 00021400 00021500 00021600 00021700 00021800 00021900 00022000 00022100 00022400 00022500 00022600 »««»»«#«»«'»««««#««««««00022700 00022800 00022900 T TERME. 00023000 PCPOINT KFB. 00023100 ; 00023200 00023300 00023400 00023500 55 BOOLEAN 8CI ; MODE := BC ; % % % % % % % % % % LOOP K := 0tLL3t BEGIN BCI := BC ; *!ODE := BCI LOOP J := 1 % BEGIN L := K TERMA : TERMC : TERMEIL COFAtLl % E := E KFAtL) % KFBtLl % KFCILl: + ( ♦ ( % KFDtLl BCI : = MODE := END : END ; MODE := TRUE ; END ; % PE THOMASONROtfS (PCPOINT IF THE NET LENGTH (RADIAL DIRECTION) IS LONGER THAN 62t THE NET MUST BE SEPARATED INTO SEVERAL EQUAL LENGTH SECTIONS. A SECTION MAY AT MOST BE 62 ROWS LONG )/W(R) END END BEGI N KFDtR) := KFCtR) LOOP WORK BACKWARD VECTORS TO GET SOLUTION X(N) DO = G(N) J := 1«1*LL5 BEGIN MODE := TRUE ; R := RTL(1»M0DE»R) BCI := REVL(1»BCI) ; MODE := BCI ; KFDIR) := KFClRl-KFBtR) *RTL(1» * X(R) = % % END ; ; ♦KFDtR-ll) ; G(R) - Q(R) « X(R+1) NEW VORTICITY VALUES STORED IN SKEWED FORM IN KFDtR) END END ; 00029600 00029700 00029800 00029900 00030000 00030100 00030200 00030300 00030400 00030500 00030600 00030700 00030800 00030900 00031000 00031100 00031200 00031300 00031400 00031500 00031600 00031700 00031800 00031900 00032000 00032100 00032200 00032300 00032400 00032500 00032600 00032700 00032800 00032900 00033000 00033100 00033200 00033300 00033400 00033500 00033600 00033700 00033800 00033900 00034000 00034100 00034200 00034300 00034400 00034500 00034600 00034700 00034800 00034900 00035000 00035100 00035200 00035300 00035400 00035500 57 MODE := TRUE ; END : PE REAL SUBROUTINE UNSKEWANDSTOREINTERMEDIATEVORTICITY (PCPOINT KED) ; BEGIN % UNSKEW RESULTS AND OVERLAY VORTICITY % WITH THE NEW VALUES CINT JfK.L ; MODE := BC ; LOOP K := 0» LL3.LL2 DO BEGIN LOOP J := 1»1«LL3 DO BEGIN L := K + J ; ZETAIL) := RTL(J-1*TRUE»KEDILJ) ; END ; END ; MODE := TRUE ; END ; NET % W t R ] : % QIR] : % GtRl : % xtR] : = KEDIR) = KFBIR] = KECIR) OR D2 OR B OR Dl = KEDIR] OR D2»W[R] 01 01 01 01 01 0( 01 01 01 0! 01 01 01 01 0! 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 0( 01 01 01 01 01 01 0< Oi Oi 01 01 01 01 01 01 01 Oi Oi PE REAL SUBROUTINE COMPUTENEWADICOEFFICIENTS (CNPOINT TERME* CNPOINT COFA* PCPOINT KEA, PCPOINT KFB» PCPOINT KECt PCPOINT KED ) ; BEGIN % DEVELOP NEW COEFFICIENTS CINT J»K»L ; PREAL TERMA* TERMC ; MODE := BC » LOOP J := 1. 1* LL1 DO 3EGIN TERMA :=-(PSIlJ+ll -PSIU-11) » TERM4I ; TERMC :=- - RTR ( 1 ♦TRUE»PSI ( J J * TERM4I ; COFAIJ] := TERMEIJ] ♦ TERMD2 ; % "A" KFAtJJ := (-TERMA+TERMB)* RTR ( 1 ♦ TRUE ♦ ZETAI J ) > ♦ (TERMEU) - TERMB2 ) » ZETAIJ] ♦(TERMA+TERMB) • RTL ( 1 ♦ TRUE* ZETA t J ] ) % ••D" ; 00049900 TRMC := ET *(PSHJ+1] ♦ PSIU-U) ; 00050000 TRMD :=( E * E * ZETAIJJ) * GT ; 00050100 % 00050300 % SAVE THE OLD VALUE OF THE STREAM FUNCTION. 00050300 % 00050400 TEST := PSHJ] ; 00050500 % 00050600 % COMPUTE NEW VALUES OF STREAM FUNCTION. 00050700 % 00050800 PSHJ) := TRMA ♦ TRMB + TRMC - TRMD ; 00050900 % 00051000 % TEST FOR CONVERGENCE 00051100 % 00051200 DIF := PSHJ] - TEST ; 00051300 IF ABS(DIF) GTR LIMIT THEN LlM := 1 5 00051400 E := E • DXS ; 00051500 END : 00051600 % o«»^o«»*«»»»'»««««««tt»#»»»*#»»**»»H»tt«#»##««»«»tt»«*««tttttttt»«*»»««'»'»*««» 00051700 % 00051800 PE REAL SUBROUTINE 00051900 STREAMFUNCTION ( PCPOINT PSI. PCPOINT ZETA ) ; 00052000 BEGIN % DETERMINE THE NEW VALUES OF THE STREAM 00052100 % FUNCTION USING THE MODIFIED SUCCESSIVE 00052200 % OVERRELAXATION (MSOR) METHOD 00052300 LABEL ITERATE ; 00052400 PREAL E ; 00052500 CINT C» PASS* LIM ; 00052600 PINT J. K. SW1. SW2 ♦ 00052700 BOOLEAN BEND. PMODE ; 00052800 % PARAMETERS REQUIRED. RLXP. DX. BC DTSI. DETA, 00052900 % Pit LL1. LL2. LIMIT.AT.BT.CT. 00053000 % ET.FT.GT.TA.DXS.PINC.LM 00053100 MODE := BC ; 00053200 NITER := ; 00053300 60 $♦ ITERATE: BEGIN 00053400 LABEL CONVERGED* PASS2 ; 00053500 % THE MODIFIED SUCCESSIVE OVERRELAXATION METHOD 00053600 % IS AN ITERATIVE METHOD WHICH REQUIRES THAT 00053700 % SUCCESSIVE VALUES OF THE STREAM FUNCTION 00053800 % CONVERGE TO THE TRUE VALUE. CONVERGENCE 00053900 % IS MEASURED BY COMPUTING THE ABSOLUTE 00054000 % DIFFERENCE BETWEEN THE VALUE OF THE STREAM 00054100 % FUNCTION ON TWO SUCCESSIVE ITERATES 00054200 % AND COMPARING IT TO AN ERROR 00054300 % LIMIT. WHEN THE ERROR IS 00054400 % LESS THAN THE LIMIT AT EACH MESH POINT, 00054500 % CONVERGENCE IS ASSUMED. LIM IS A CU 00054600 % VARIABLE WHICH IS SET TO ZERO AT THE 00054700 % BEGINNING OF EACH ITERATION. IF ANY MESH 00054800 % POINT FAILS THE CONVERGENCE TEST, LIM 00054900 % IS THEN SET EQUAL TO ONE AND THUS 00055000 % INDICATES TO THE PROGRAM THAT ANOTHER 00055100 % ITERATION IS NEEDED. 00055200 LIM := o ; 00055300 % DEVELOP THE TERMS IN THE FINITE DIFFERENCE 00055400 % FORMULA FOR FINDING AN IMPROVED VALUE 00055500 % OF THE STREAM FUNCTION. 00055600 DISPLAY 0,11 00055700 PASS := 1 ; 00055800 SW1 := CEVEls i ; 00055900 SW2 := CODD ; 00056000 BEND := BOD ♦ 00056100 NITER := NITER + 1 ; 00056200 PASS?: FOR ALL SW1 = 1 DO E := PI » DXS ; 00056300 FOR ALL SW1 = DO E := PI * DX ; 00056400 LOOP C := li 2» LM DO 00056500 BEGIN 00056600 J := C + swi ; 00056700 K := C + sw2 ; 00056800 PMODE : = mode ; 00056900 NEWPSI (J,K,LIM,E»PMODE) ; 00057000 MODE := : bc ; 00057100 END ; 00057200 % IF THERE ARE AN ODD NUMBER OF ROWS IN THE 00057300 % PSI NET, THE LAST ROW MUST BE UPDATED BY 00057400 % THE FOLLOWING PROCEDURE WHICH OPERATES 00057500 % ON ONLY THE ODD COLUMNS IN PASS ONE AND 00057600 IF EVEN = THEN 00057800 % THE EVEN COLUMNS IN PASS TWO. 00057700 BEGIN 00057900 J := LLl 5 00058000 K := LLl ; 00058100 PMODE : = BEND ; 00058200 NEWPSI (J»K,LIM,E»PMODE) ; 00058300 MODE := s bc ; 00058400 END i 00058500 % AFTER PASS ONE THE CONTROL VARIABLES ARE 00058600 % ADJUSTED FOR PASS TWO. 00058700 IF PASS = 1 THEN 00058800 BEGIN 00058900 PASS := > 2 ; 00059000 SW1 : = codd ; 00059100 SW2 : = CEVEN 1 00059200 BEND := : BEV ; 00059300 6l GO TO PASS2 ; 00059400 END ; 00059500 IF LIM = THEN GO TO CONVERGED ; 00059600 IF NITER = 6 THEN GO TO CONVERGED ♦ 00059700 GO TO ITERATE ; 00059800 CONVERGED: MODE := TRUE ♦ 00059900 END ; 00060000 END ; 00060100 % 00060300 PE REAL SUBROUTINE 00060400 SURFACEVORTICITY (PCPOINT PSI, PCPOINT ZETA) ; 00060500 % COMPUTE THE SURFACE VORTICIY 00060600 % USING A TAYLOR SERIES APPROXIMATION 00060700 BEGIN 00060300 CREAL PRODUCT, PREFIX ; 00060900 % PARAMETERS REQUIRED, DTSI, PI 00061000 MODE := BCS ; 00061100 PRODUCT := PI • DTSI ; 00061?00 PREFIX := 1.0 / ( 2.0 » PRODUCT * PRODUCT ) ; 00061300 ZETA101 := PREFIX * ( 8.0 * PSItl] - PSK2) ) ; 0006140C MODE := TRUE ; 00061500 END ; 00061600 % «■»<*«•«•»•» tt«»tt«tt»#«**»#*»»#**»» ^•»«-»tt»»««'<»»»#^<*'tt<>tt«tttt'innnni- ««««»»««••»««•«<)•<» 00061700 % 00061800 PE REAL SUBROUTINE 00061900 INITIALIZE ( PCPOINT PSI, PCPOINT ZETA ) ; 00062000 BEGIN 00062100 CREAL EX, IX ; 00062200 CINT JD, JDE ; 00062300 % PARAMETERS REQUIRED, DX, IDX, LLl, S 00062400 % 00062500 BEGIN % USE THE POTENTIAL FLOW SOLUTION AS THE 00062600 * INITIAL STREAM FUNTION VALUES AND 00062700 % ZERO THE VORTICITY NET 00062800 MODE := BCS ; 00062900 PSHOJ := 0.0 I 00063000 ZETAtO] := 0.0 ? 00063100 EX := DX ; 00063200 IX := IDX ; 00063300 LOOP JD := 1,1, LLl DO 00063400 BEGIN 00063500 MODE := BCS ; 00063600 ZETAIJD1 := 0.0 ; 00063700 PSKJD] := 0.0 ; 00063800 MODE := BC ; 00063900 PSItJD] := (EX - IX)»S ; 00064000 EX:= EX * DX ; 00064100 IX:= IX * IDX ; 00064200 JDE := JD ; 00064300 END ; 00064400 JD := JDE + 1 ; 00064500 PSItJDl := EX « S ; 00064600 ZETACJD1 := 0.0 ; 00064700 MODE := TRUE ; 00064800 END 5 00064900 PRINT1 (PSI,3) ; 00065000 STREAMFUNCTION ( PSI, ZETA ) ; 00065100 SURFACEVORTICITY ( PSI, ZETA ) ; 00065200 END; 00065300 62 % 00 PE REAL SUBROUTINE PRINTPARAMETERS » BEGIN MODE := TRUE ; SIMWRITE (LINE. "REYNOLDS NUMBER'S NREY. ••NUMBER OF CIRCUMFERENTIA •'NUMBER OF RADIAL MESH PO ••CIRCUMFERENTIAL INCREMEN "RADIAL INCREMENT SIZE"»D "MAXIMUM CIRCUMFERENTIAL "MAXIMUM RADIAL COORDINAT "DIMENTIONLESS TIME INCRE "MAXIMUM DIMENTIONLESS TI TLMT. "NUMBER OF TIME STE PINC. "RELAXATION PARAMETE "ERROR LIMIT FOR ITERATIV LIMIT. "BOOLEAN MASKS FOR BCS.BCT) ; END ? % PE REAL SUBROUTINE OUTPUTRESULTS (CINT OUT ST. CREAL T) ; BEGIN MODE := TRUE ; IF ST EQL PINC THEN BEGIN SlMWRITE(LINE."DIMENSIONLESS TIME".T ) ; ORINT1 (ZETA.l) ; ORINTl (PSI.2) ; SIMWRITE (LINE. "NUMBER OF ITERATIONS = " . NITER M ST := ; DISPLAY OtlS END ; ST := ST + l ; END ; % L MESH POINTS". NETA INTS".NTSI» T SIZE".DETA, TSI. COORDINATE".ETAMAX. E'STSIMAX. MENT"»DT. ME DESIRED". PS BETWEEN PRINT-OUT R"»RLXP. E CONVERGENCE". FIELD WIDTH'S »e«tt#tt«a«tttt«tt«e«tttttttt $• MAIN BEGIN LABEL REPEAT ; CREAL T ; CINT ST ; PE REAL SUBROUTINE * ACCUMULATED DIMENSIONLESS TIME % PRINTOUT PARAMETER MODE := TRUE ; T := 0.0 ; ST := 1 ; SYSTEMPARAMETERS ; COMPUTESYSTEMCONSTANTS ; PRINTPARAMETERS ; INITIALIZE ( PSI. ZETA ) ; OUTPUTRESULTS (ST.T) ; % INITIALIZES THE STREAM FUNCTION AND VORTICI % NETS. RELAXES THE PSI NET. AND COMPUTES % INITIAL SURFACE VORTICITY REPEAT:T := T + DT S 00 00 00 00 00 ♦ 00 00 00 00 00 00 00 00 S".00 00 00 00 00 00 »00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 •»»»oo 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 TY 00 AN 00 00 00 065400 065500 065600 065700 065800 065900 066000 066100 066200 066300 066400 066500 066600 066700 066800 066900 067000 067100 067200 067300 067400 067500 067600 067700 067800 067900 068000 068100 068200 068300 068400 068500 068600 068700 068800 068900 069000 069100 069200 069300 069400 069500 069600 069700 069800 069900 070000 070100 070200 070300 070400 070500 070600 070700 070800 070900 071000 071100 071200 071300 63 % INCREMENT TIME 00071400 VORTICITY ( PSI. ZETA ) ; 00071500 % COMPUTE NEW VALUES OF VORTICITY BY ADI 00071600 STREAMFUNCTION ( PSI* ZETA ) ; 00071700 % COMPUTE NEW VALUES OF THE STREAM FUNCTION 00071800 % BY THE MODIFIED SUCCESSIVE OVERRELAXAT ION 00071900 % METHOD 00072000 IF ST NEQ PINC THEN 00072100 SIMWRITE (LINE* "NUMBER OF ITERATIONS = "» NITER M 00072200 SURFACEVORTICITY ( PSI. ZETA ) ; 00072300 % ESTIMATES SURFACE VORTICITY BY REFERENCE TO 00072400 % THE STREAM FUNCTION VALUES NEAR THE SURFACE00072500 OUTPUTRESULTS (ST»T) ; 00072600 IF T LSS TLMT THEN GO TO REPEAT ; 00072700 END ; 00072800 % *-tHHHHHHHHHHHHHHHHHHH» «««««««««#«««##»««tt»»»«tf»»»«4ttttt«»Utt«««4ftOtt«tt««ft()Q0 72900 % tt»»*tt»«»»»»»»»»»*«»*#«»*«#*»«-#»»»*«*»»»#«»»»«»»«»«#»#'»»»*'tt#»#«»»Hi-»'»»oOO73OO0 % 00073100 MAIN ; 00073200 FINI: END . 00073300 6k APPENDIX B 4ZIP JLI«RARY * CICHY/GLYOIRECT % 5/24/72 BEGIN CREAL Pit OETA. OTSI» DT« DX» IDX. NREYf KE» TLMT* ETAMAX. TSIMAX* TERMM» TERMB* TERMR2. TER*D.TERMD2.DTI2»HXtHY t CINT META. NTSIt LL 1 * LL2»LL3»LL4.LL5»LL6» NITER. PINO CN»CMtPP»I.J«CK.CS*OFr.PPo«K t PINT CEVEN» CODQJ PINT VECTOR INDEX! 1) I PREAL ETA. M. $. Xt E » PREAL VECTOR PSK 20). ZETA( 20 J . AC 1 1 .SSt 1 J . ISC 1 J BOOLEAN RCS1. BC11. BCS. SAVEMODE. 8C AD.BCT ; FILE LINE (10 ROWS ) I FILE UOISK1 ="A / OIS » LABEL FINI I ; % % % % % % *♦ PE REAL SUBROUTINE EXP AS RGA (PE REAL ARG AS RGA) I % EXPONENTIAL BEGIN CALL EXP64 I END » 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 • *>«»#iHHKHHKHH»00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 • •••«4»««*«»»«00 00 00 00 00 00 00 00 •••••••••••»«00 00 00 00 00 00 00 PE REAL SUHROUTINE SIN AS RGA (PE REAL ARG AS RGA ) I % SINE BEGIN CALL SIN64 ; END « PE REAL SUBROUTINE COS AS RGA (PE REAL ARG AS RGA); % COSINE BEGIN CALL COS64 I END i PE REAL SUBROUTINE LN AS RGA ( PE REAL ARG AS RGA ) ; % NATURAL LOGARITHM BEGIN i° CALL LN64 ; END I % ««««««««»«»«»«»««»««»»««»««»«t»*«0»«*»*«»**«»»»*t ; 00016400 3CS := BCS OR BC11 ; 00016500 END ; 00016600 % THE BOOLEAN VARIABLE BC MUST REFLECT THE 00016700 % WIDTH OF THE NET. A "1" SHOULD BE PLACED 00016300 % IN EACH BIT POSITION WHICH CORRESPONDS 00016900 * TO A PE CONTAINING ONE OF THE INTERNAL MESH00017000 % POINTS. ALL OTHER BITS MUST BE ZERO. 00017100 BEGIN 00017200 CINT j; 00017300 BCT := FALSE ; 00017400 LOOP J := 1»1«LL3 DO 00017500 BEGIN 00017600 BCT : = BCT OR BCS1 ; 00017700 BCT : = REVR(1,BCT) ; 00017800 end; 00017900 end; 00018000 TERM4I := 0.25 /(DTSI » DETA) ; 00018100 67 TERMB TERMB2 TERMD TERMD2 DTI2 = 2.0 /(NREY • DETA * DETA) ; - 2.0 » TERMB ; = 2.0 /(NREY • DTSI • DTSI) ; = 2.0 • TERMD ; = 2.0/OT ; % 00018200 00018300 00018400 00018500 00018600 00018700 MODE := TRUE ; MODE := (PEN LEQ (NETA - 1) / 2 ) ; ETA := PEN ; MODE := TRUE ; MODE := (PEN GTR (NETA - 1) / 2 ) ; ETA := (NETA - 1 ) - PEN ; MODE := BC ; M := ETA t PI * DETA ; S := SIN(M) ; MODE := TRUE ; END ; % PE REAL SUBROUTINE COMPUTEADICOEFFICIENTSANDSKEW (CREAL E» CNPOI CNPOINT COFA» PCPOINT KFA PCPOINT KFC» PCPOINT KFD BEGIN % SKEW COEFFICIENTS CINT K»J*L ; PREAL TERMA* TERMC ; BOOLEAN 8CI ; MODE := BC ; % IF THE NET LENGTH (RADI NT TERME. » PCPOINT KFB. ) ; % LONGER THAN 62* THE % INTO SEVERAL EQUAL L % SECTION MAY AT MOST % (LL3=62> AND 62 WIDE % PROGRAM TREATS ONE S % ALL SECTIONS HAVE BE % NEW VORTICITY VALUES % PROGRAM WHEN THE LAS % PROCESSED. LOOP K := 0.LL3. LL2 DO BEGIN 3CI := BC ; viODE := BCI ; LOOP J := It 1» LL3 DO % THERE ARE LL3 ROUTES IN % SECTION TO BEGIN L := K ♦ J ; TERMA :=-(PSHL + l) -PSKL-ll) TERMC :=-(RTL(l»TRUE»PSI[L))-R • TERM4I ; TERMEIL1 := E»E»DTI2 ; COFAILJ := TERMECL] + TERMB2 ; % " A •' AL DIRECTION) IS NET MUST BE SEPARATED ENGTH SECTIONS. A BE 62 ROWS LONG (LL4=62) . THE ECTION AT A TIME UNTIL EN SOLVED FOR THE . LL2 TELLS THE T SECTION HAS BEEN EACH SKEW THE ARRAY. * TERM4I ; TR(1,TPUE»PSI(LJ)> 00018800 00019100 00019200 00019300 00019400 00019500 00019600 00019700 00019800 00019900 00020000 00020100 00020200 00020300 00020400 00020500 00020600 00020700 00020800 00020900 00021000 00021100 00021200 00021300 00021400 00021500 00021600 00021700 00021800 00021900 00022000 00022100 00022200 00022300 00022400 00022500 00022600 00022700 00022800 00022900 68 ♦TERMA - TERMB ) ; ♦-TERMA - TERMB ) i E := E * DX ; KFA(L) := RTR(J-1» % •• C •• KFBtLl := RTR(J-1, % » B " KFC[L):= RTR(J-1» ♦ (TERMC+TERMD) * ZETAtL-1] +(TERMEIL1-TERMD2) * ZETACU +<-TERMC+TERMD) « ZETAIL+11); % '• D " KFDIL) := C0FAIL1 ; BCI := REVR(1»BCI) ; MODE := BCI ; END ; END ; MODE := TRUE ; END ; % PE REAL SUBROUTINE THOMASONROWS (PCPOINT KFA, PCPOINT KFB» PCPOINT KFC» PCPOINT KFD ) BEGIN % APPLY THE THOMAS METHOD TO EACH ROW IN THE % SECTION SIMULTANEOUSLY THUS INVERTING % LL3 MATICIES AT A TIME. % % REFERENCE TEXT: LAPIDUS, L. . DIGITAL % COMPUTATION FOR CHEMICAL ENGINEERS* % MC GRAW-HILL* 1962» P. 254. CINT J»K ; MODE := TRUE ; LOOP K := 0»LL3» LL2 BEGIN PINT R ; BOOLEAN BCI ; BEGIN CINT N ; DO BCI := BCT; MODE := TRUE ; R := PEN + K ; MODE := BCI ; % WORKING FORWARD COMPUTE % INTERMEDIATE VALUES % BC CARRIES THE BIT PATTERN % REFLECTING THE WIDTH OF THE % NET. THIS INFORMATION IS % PASSED TO BCI FOR USE IN % ENABLING THE PROPER % PROCESSING ELEMENTS WHILE % OPERATING ON THE SKEWED % APRAY. % R IS LIKE A POINTER WHICH % POINTS TO THE PROPER % VECTOR ELEMENT IN EACH PE % W(l) = A(l) KFB[R] := KFBfR) / KFDtR] * % Q(l) = B(l) / W(l) KFCtRl := KFCIR] / KFDtR] ; % 6(1) = D(l) / W(l) LOOP J := 2»1»LL4 DO BEGIN MODE := TRUE ; % SHIFT INDEXING AND ENABLING 00023000 00023100 00023200 00023300 00023400 00023500 00023600 00023700 00023800 00023900 00024000 00024100 00024200 00024300 00024400 00024500 »00024600 00024700 00024800 ;00024900 00025000 00025100 00025200 00025300 00025400 00025500 00025600 00025700 00025800 00025900 00026000 00026100 00026200 00026300 00026400 00026500 00026600 00026700 00026800 00026900 00027000 00027100 00027200 00027300 00027400 00027500 00027600 00027700 00027800 00027900 00028000 00028100 00028200 00028300 00028400 00028500 00028600 00028700 00028800 00028900 69 % PARAMETERS r := RTR(1.M0DE»R) ; BCi:= REVR(1»BCI) ; MODE := 8CI ; KFD(R]*.= KFDtRl - KFACRJ *RTR(1» »KFB(R+11> ; % W(R) = A(R) - C(R» • Q(R-1 KFBtRJ := KFBCR1 / KFD[R] ; * O(R) = B(R) / W(R) KFCtR]:=(KFC(Rl- KFAtR) * RTR(1» ,KFC[R+1J) ) / KFDtRl ; % G(R)=( D(R)-C(R)*G(R-1) )/ END * END ; BEGIN KFDtRl := KFCIR] ; % WORK BACKWARD TO GET SOLUT % VECTORS % X(N) = G(N) DO ION LOOP J := 1»1»LL5 BEGIN MODE := TRUE ; R := RTL(l*MODE»R) ; BCI := REVL(l.BCI) ; MODE := BCI ; KFD[R] := KFCIRJ-KFBIR] »RTL<1» .KFDIR-11) ; % X(R) = G(R) - Q(R) * X(R+1 % NEW VORTICITY VALUES STORE % IN SKEWED FORM IN KFD[R END ; END ; END ; MODE := TRUE ; END ; * PE REAL SUBROUTINE UNSKEWANDSTOREINTERMEDIATEVORTICITY (PCPOINT KFD) ; BEGIN % UNSKEW RESULTS AND OVERLAY VORTICITY NET % WITH THE NEW VALUES CINT J»K»L ; MODE := BC ; LOOP K := 0, LL3»LL2 DO 9EGIN LOOP J := ltl.LL3 DO BEGIN L := K ♦ J ; ZETAIL] := RTL(J-l»TRUEtKFDCL]) ; END ; END ; MODE := TRUE • END ; % W[R) : •= KFDtRl OR D2 % Q[R] :• := KFBCR] OR B % GIR] :, = KFCIR] OR Dl % X[R) :: = KFDtR) OR D2tW[R] PE REAL SUBROUTINE ) D ] 0029000 0029100 0029200 0029300 0029^00 0029500 0029600 0029700 0029800 0029900 0030000 0030100 0030200 0030300 0030400 0030500 0030600 0030700 0030800 0030900 0031000 0031100 0031200 0031300 0031400 0031500 0031600 0031700 0031800 0031900 0032000 0032100 0032200 0032300 0032400 0032500 0032600 0032700 0032800 0032900 0033000 0033100 0033200 0033300 0033400 0033500 0033600 0033700 0033800 0033900 0034000 0034100 0034200 0034300 0034400 0034500 0034600 0034700 0034800 0034900 70 COMPUTENEWADICOEFFICIENTS (CNPOINT TERME* CNPOINT COFA» PCPOINT KFA, PCPOINT KFB» PCPOINT KFC» PCPOINT KFD ) ; BEGIN % DEVELOP NEW COEFFICIENTS CINT J»K»L ; PREAL TERMA* TERMC ; MODE := BC ; LOOP J := 1» 1. LL1 DO BEGIN TERMA :=-I / NX ; H := HX / HY ; KN := NX.t 16:411 ; LOOP I := 0*1*KN DO BEGIN ^ODE := TRUE ; MODE := (L LEQ LL6 / 2 LX := l ; CO := COS(LX * X) ; MODE := TRUE ; MODE := (L GTR LL6 / 2 LX := LL6 - L ; CO := -COS(LX * X) ; MODE := TRUE ; Sill := SIN(LX»X) ; 00045400 00045500 00045600 *»*00045700 00045800 00045900 00046000 00046100 00046200 00046300 00046400 00046500 00046600 00046700 00046800 00046900 00047000 00047100 00047200 00047300 00047^*00 00047500 »*»00047600 00047700 OOO47RO0 00047900 00048000 00048100 00048200 00048300 00048400 00048500 00048600 00048700 00048800 X. 00048900 S 00049000 00049100 00049200 00049500 00049600 00049700 00049800 00049900 00050000 00050100 00050200 00050300 73 ISII) := sen +stn ; IF ISC I J NEQ THEN BEGIN ISI 1 1 := -1.0 / IStll ; END; All) := -2.0 * (H-H * CO ♦ 1 ) ; INDEXU1 := TURN(PP-l.L) ; L := L + 64 ; END ; END ; % SUBROUTINE FOURIER(PC°OINT P»CINT K.CINT CNtCINT LN,BOOLEAN FMODE. PCPOINT DS»PCPOINT IS) ; BEGIN CINT Jl ; BOOLEAN SMODE ; CINT ItL»J»IL»N2»ILJ»CJ.KM,N3»IT,IP,Nll.I2.Il»IM.IQ,KP.KMl» PINT »II PREAL 0DDlt0DD2»0DD3; CREAL RR»RR1»C»S ; J :=(CN-1)»K; L:=J.[ 16:471; ppp:=i; PI :=o; IF ( PEN EQL OFF+1 ) THEN PI : = -1 ; MODE := FMODE ; KM1 := K-l ; kp := k + K ; N3 := CN - 1 ; N2 := N3.[ 16:471; LOOP ILJ := Otl.KMl DO 3EGIN km := k + ILJ ; LOOP I := K»K»L-K DO BEGIN IL := I - ILJ ; ODD1 := P[ I+ILJ+PIi; PtPI+I+ILJ] := ODD1 + PCPI+J-IL1 ; PfPI+J-ILl := ODD1 - PIPI+J-IL1 ; end; d[PI+L+ILJ1 := PtPI+L+ILJl + PtPI+L+ILJl ; IL := K - ILJ ; OCPI+ILJ1 := - PCPI+J-IL1 - PCPI+J-IL1 ; 0DD3 := P[PI+KM]; PtPI+KM) := -0DD3 - ODD3 ; IL := ; LOOP I := KP.KP»L-KP DO BEGIN IL := IL ♦ 2 ; CJ := I ♦ ILJ ; ODD1 := PCPI+I+KM1 - 0DD3 ; GL := IL .116:421 ; SMODE := MODE ; MODE := TRUE ; RR := GRAB0NEP I := KP»KP,J-KP DO 00058300 BEGIN 00058400 IL := IL + 1 ; 00058500 GL := IL .[16:42] ; 00058600 0DD3 := INDEXIGL1 ; 00058700 CJ := GRAB0NE(0DD3 *IL) ; 00058800 IF IL LSS CJ THEN 00058900 BEGIN 00059000 CJ := CJ * KP ; 00059100 LOOP ILJ := 0»1»KM DO 00059200 BEGIN 00059300 ODD1 := PI I+ILJI ; 00059400 P[ I+ILJJ := P [ CJ+ILJ1 ; 00059500 P[ CJ+ILJ) := ODD1 ; 00059600 END ; 00059700 END ; 00059800 END ; 00059900 IT := kp + kp ; 00060000 IP := N3 ; 00060100 LOOP I := 2»1»LN -1 DO 00060200 3EGIN 00060300 IP:=IP.[ 16:47] ; 00060400 Nil := o ; 00060500 12 := l ; 00060600 CJ := IT. (16:47) ; 00060700 IT := IT + IT ; 00060800 LOOP IL := 0»CJ»CJ DO 00060900 BEGIN 00061000 LOOP KM := IL»KP*IL+CJ-KP DO 00061100 BEGIN 00061200 GL := (N2-N11) .116:42) ; 00061300 C := GRABONE(DSIGL)»N2-Nll) * 12 ; 00061400 GL := Nil .[ 16:42) ; 00061500 S := GRABONE(DS(GL)*Nll) ; 00061600 LOOP Jl := OtIT.J-IT DO 00061700 BEGIN 00061800 MODE := FMODE ; 00061900 LOOP ILJ := 0.1»KM1 DO 00062000 BEGIN 00062100 IM := Jl + KM + ILJ • 1 00062200 IQ := IM + CJ ♦ CJ \ 00062300 ODD1 := Pf PI+IQ)»C-P[PI+IQ+K] » S ; 00062400 75 END ; := Nil 0DD2 := PtPI + IQ+K)»OP[PI + IQ] • S PIPI+IQ1 := PtPI+IMl - 0DD1 I PCPI+IQ+KJ := PCPI+IM+KJ - 0DD2 ; PIPI+IM1 := PIPI+IM] + ODD1 ; PtPI+IM+K] := PtPI+IM+K) ♦ 00D2 ; MODE := TRUE ; END ; + IP ? Nil END ; 12 := -1 ; IP I* - IP I end ; end ; I := o ; LOOP II := 1»1*N2 DO BEGIN I := I + K ; 3L := II .116:42) ; PR := GRABONE(IS(GLJ»Il) ; MODE := FMODE ; LOOP ILJ := OtltKMl DO BEGIN ODD1 := PtPI+I+ILJJ - PtPI+J-I+ILJ] ? 0DD2 := ( PCPI+I+ILJ] + PCPI+J-I+ILJ ) ) * RR ; PIPI+I+ILJ1 := 0DD2 + ODD1 ; PIPI+J-I+ILJ1 := 0DD2 - ODD1 ; MODE := TRUE ; END ; END ? END ; % POISSONDIRECT SUBROUTINE (PCPOINT P»PCPOINT CtPCPOINT A»PCPOINT S.CINT CINT CM»CREAL HXtCREAL HY, PCPOINT IS* CINT LN ) ; FUNCT ION %%m%%%%%%%%%%%%%%%%*>%% : *>%%%%<*>%%%% CN POISSON 1 EQUATION USING A DIRECT % THIS SUBROUTINE SOLVES METHOD SIMILAR TO HOCKNEY'S PROGRAM [1J. CfI»Jl = HX»(P[I*l»J)*P[I-lf J)-2*PII»J1 ♦ HY»(PtI»J+ll+PII»J-ll-2»Pt I»J1) . THE CHARGE* C» AND THE BOUNDARY POINTS OF THE POTENTIAL* P. ARE USED TO CALCULATE THE INTERIOR POINTS OF POTENTIAL. PARAMETERS %%%%%%%%%%%%%%%%%%%%i>%%%%%%%i>%%%%% % P IS THE POTENTIAL. C IS THE CHARGE. A CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRICES SOLVED IN THE SECTION OF THIS ROUTINE CALLED CRED. S CONTAINS THE VALUES OF THE SIN WHICH ARE NEEDED BY THE FOURIER TRANSFORM PART OF THE ALGORITHM. IS IS ALSO NEEDED BY FOURIER. INDEX IS GLOBAL FROM THE MAIN FOURIER. A* S» IS* AND INDEX P HAS CN ROWS AND CM COLUMNS. PROGRAM AND IS ARE CALCULATED HX(HY) IS THE VERTICAL) COEFFICENT OF THE LOG BASE 2 OF CN-1. POISSON'S EQUATION. NEEDED BY BY SETPT1. HORIZONTAL ( LN EQUALS ; o o o o o o o o o o o o o o o o o 0062500 0062600 0062700 0062800 0062900 0063000 0063100 0063200 0063300 0063400 0063500 0063600 0063700 0063800 0063900 0064000 0064100 0064200 0064300 0064400 0064500 0064600 0064700 0064800 0064900 0065000 0065100 0065200 0065300 0065400 0065500 0065600 0065700 0065800 0065900 0066000 0066100 0066200 0066300 0066400 0066500 0066600 0066700 0066800 0066900 0067000 0067100 0067200 0067300 0067400 0067500 0067600 0067700 0067800 0067900 0068000 0068100 0068200 0068300 0068400 76 % 00068500 %%%%%%%%%%%%%%*>% METHOD %%%%**%%%%%%%%%%%%*<**%*%**%%%%%%%%%%*% 00068600 % 00068700 % A FOURIER ANALYSIS IS DONE TO EACH COLUMN BY SUBROUTINE 00068800 % FOURIER. THIS IS FOLLOWED BY SOLVING THE TRIDIAGONAL 00068900 % MATRICES. ONE PER ROW. FINALLY A FOURIER SYNTHESIS IS 00069000 % PERFORMED ON EACH ROW TO GIVE ANSWER. ■ FOR FURTHER DETIAL 00069100 % SEE [21. 00069200 % 00069300 %%%%%%%%%%%%%%*%%%%%% REFERENCES %%%%%%%*%%**%%*%%*%%%%%%%%%%% 00069400 BEGIN 00069500 CINT K»CM2»CN2«CMIK»CNIK»CI»CJ« CL» IL.L1 ♦L2»L3»L4 »KM1»KP»KC; 00069600 PINT d i*pj; 00069700 CREAL POTFAC; BOOLEAN FMODE ; 00069800 PCPOINT STORE ; 00069900 PREAL PR» E» Q* T ; 00070000 MODE := TRUE ; 00070100 CMIK := CM .[ 16:42] ; 00070200 CI := ; 00070300 IF CM. [58:61 GTR THEN CI := 1 • 00070400 CN2 := CN - 2 ; POTFAC := 0.125 / < (CN- ■1)»HY) ; 00070500 CM2 := CM - 2 ; 00070600 K := CMIK + CI ; 00070700 kmi := k - l ; 00070800 kp := k + k ; 00070900 STORE := GETPEB(KP+K) ; 00071000 KC := K * CN2 ; PI := o ; 00071100 CI := CM. [58:6] ; 00071200 OFF := CI - 2 ; 00071300 IF CI = THEN OFF := 62 ; 00071400 MODE := ( PEN GTR OFF ) ; 00071500 PI := -1 ; 00071600 MODE := TRUE ; 00071700 CJ := KC + 1 ; 00071800 LOOP CI := 0»1»KM1 DO 00071900 BEGIN 00072000 STOREtCIl := CtCI+K] ; 00072100 STOREICI+K] := CICI+KC1 ; 00072P00 STORE (CI+KP] := P[CI] ; 00072300 77 END ; LOOP CI := 1»1»KM1 DO BEGIN CIPI+CJ1 := CIPI+CJ] - PCCJ+PI+K1 * HX ; CJ := CJ +1 ; CIPI+K+CI) := CIPI+K+CI1 - P[PI+CI1 » hx ; end; if pen gtr and pen lss off + 1 then 3EGIN CtK] := CCK] -PtO] * hx ; CtKC] := CIKCJ - P[KC+K] • HX ; FMODE := MODE ; END ; pj := ; PI := ; MODE := ( PEN EQL OFF ) ; PI := kmi ; MODE := TRUE ; MODE := ( PEN GTR OFF-1 ) ; PJ := -1 ; MODE := TRUE ; IF OFF = 1 THEN BEGIN FMODE :=(PEN GTR 0); END; LOOP CI := K»K»KC DO BEGIN 3R := 0.0 ; MODE := ( PEN EQL OFF ) ; »R := RTL(1»»PCCI+KM1 )) ; MODE := TRUE ; *!ODE := PEN EQL 1 ; PR := RTR(1»»P£CI)) + PR ; viODE := TRUE ; MODE := FMODE ; = > [CI+PI1 : = (C[CI + PI) -PR*HY) * POTFAC ; MODE := TRUE I LOOP Ll:= 1»1»KM1 DO BEGIN PtLl+CI+PJ) := CtLl+CI+PJ] » POTFAC ; end; end ; fourier ; AtCIl := RTR(lt»AtCIl) ; E := 1.0 / ( AtCIl - RTRUttE) ) PI := RTR(1»«PI) ; Q := (PIPI+L1+L4] - RTRUttQ)) • L4 := ; END ; L4 := -1 ; IL := 1 ; end; CL := 64 ; E := PtPI+KMl) - AtCIl * Q ; PIPI+KMl l := q ; IL := 4 ; L4 := o ; LOOP LI := 0»lfKMl DO BEGIN IF L1=KM1 THEN CL := CJ ; L3 := KM1 - LI ; LOOP L2 := IL«1»CL DO BEGIN MODE := REVLCltMODE) ; PI := RTL(l»tPI) ; A(CI 1 := RTL(1».A[CI1) ; IL END = PI Q := RTLU»»Q) ; E := RTL(UtE) ; T := PCPI+L3+L4] Q := E ; E := T ; PIPI+L3+L4] := Q L4 := ; END ; = l ; L4 := 1 ; ♦ + 64 * k ; - Q - AtCI 1 * E I TRUE ;<* UNSKEW POTENTIAL KP»K,KC DO i ; CJtl»CJ+KMl DO »I END MODE ;= CI := ; LOOP CJ := BEGIN CI := CI ♦ LOOP IL := BEGIN PI ILl := RTL(CI»»P( IL1) END ; END ; FOURIER (P» K»CN»LN»FMODE»S» LOOP CI := 0»1»KM1 DO BEGIN CtCI+K] := STORECCI 1 ; CtCI+KC) := STOREtCI+K] ; IS 00078400 00078500 00078600 00078700 00078800 00078900 00079000 00079100 00079200 00079300 00079400 00079500 00079600 00079700 00079800 00079900 00080000 00080100 00080200 00080300 00080400 00080500 00080600 00080700 00080800 00080900 00081000 00081100 00081200 00081300 00081400 00081500 00081600 00081700 00081800 00081900 00082000 00082100 00082200 00082300 0OO824O0 00082500 00082600 00082700 00082800 00082900 00083000 00083100 00083200 00083300 00083400 00083500 00083600 00083700 00083800 00083900 00084000 00084100 00084200 00084300 79 p(CI] := STOREtCI+KP] ; end ; 00 00 end ; 00 % PE REAL SUBROUTINE SURFACEVORTICITY (PCPOINT PSI, PCPOINT ZETA) I % COMPUTE THE SURFACE VORTICIY % USING A TAYLOR SERIES APPROXIMATION BEGIN CREAL PRODUCT* PREFIX ? % PARAMETERS REQUIRED* DTSI, PI MODE := BCS ; PRODUCT := PI » DTSI ; PREFIX := 1.0 / ( 2.0 • PRODUCT * PRODUCT ) ; ZETAfOJ := PREFIX • ( 8.0 * PSIUJ - PSII21 ) ; MODE := TRUE ; END ; 00 00 00 00 00 00 00 00 00 00 00 00 00 00 % »w#»#»««#»»»»»*»«*»*»»»«#«»»«#»'»»»##»»«*tt»»#*»»**«*«»«'»»»*##'»#*«'»-t»<»»00 % PE REAL SUBROUTINE INITIALIZE ( PCPOINT PSI, PCPOINT ZETA ) ; BEGIN CREAL EX, IX ; CINT JD, JDE ; % PARAMETERS REQUIRED, DX, IDX, LL1, S % BEGIN % USE THE POTENTIAL FLOW SOLUTION AS THE % INITIAL STREAM FUNTION VALUES AND % ZERO THE VORTICITY NET MODE := BCS ; pshoi := o.o ; ZETAI01 := 0.0 ; EX := DX ; IX := idx ; LOOP JD := 1,1, LL1 DO BEGIN MODE := BCS t ZETAtJDl := 0.0 ; PSHJDJ := 0.0 ; MODE := BC ; PSIIJD1 := (EX - IX)»S ; EX:= EX * DX ; IX:= IX * IDX ; JDE := JD ; END ; JD := JDE ♦ 1 I oSHJDl := EX * S ; ZETACJD1 := 0.0 ; MODE := TRUE ; END ; print1 (psi,3) ; poiss0ndirect(psi,zeta,a,ss,ntsi,neta,hx,hy,is,pp) ; surfacevorticity ( psi, zeta ) ; end; % 00 PE REAL SUBROUTINE 00 PRINTPARAMETERS ; 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 BEGIN 00 084400 084500 084600 084700 084800 084900 085000 085100 085200 085300 085400 085500 085600 085700 085800 085900 086000 086100 086200 086300 086400 086500 086600 086700 086800 086900 087000 087100 087200 087300 087400 087500 087600 087700 08780C 087900 088000 088100 088200 088300 088400 088500 088600 088700 088800 088900 089000 089100 089200 089300 089400 089500 089600 089700 089800 089900 090000 090100 090200 090300 % s» MODE := TRUE ; SIMWRITE (LINE. "REYNOLDS NUMBER'S NREY, '•NUMBER OF CIRCUMFERENTIAL MESH POINTS". NET "NUMBER OF RADIAL MESH POINTS" ,NTSI , "CIRCUMFERENTIAL INCREMENT SIZE'SDETA, "RADIAL INCREMENT SIZE"»DTSI» "MAXIMUM CIRCUMFERENTIAL COORDINATE" tETAMAX ♦ "MAXIMUM RADIAL COORDINATE" ♦ TSIMAX, "DIMENTIONLESS TIME INCREMENT" ,OT, "MAXIMUM DIMENTIONLESS TIME DESIRED"* TLMT, "NUMBER OF TIME STEPS BETWEEN PRINT-OU PINC, "BOOLEAN MASKS FOR FIELD WIDTH", BCS»BCT) ; END * a********************************************* ****************** PE REAL SUBROUTINE OUTPUTRESULTS (CINT OUT ST, CREAL T) ; BEGIN MODE := TRUE ; IF ST EQL PINC THEN BEGIN SIMWRITE(LINE»"DIMENSIONLESS TIME"»T ) ; 3RINT1 ; nnrT ^ % INITIALIZES THE STREAM FUNCTION AND VORTIC % NETS, RELAXES THE PSI NET, AND COMPUTES % INITIAL SURFACE VORTICITY REPEATtT := T + DT ; % INCREMENT TIME VORTICITY ( PSI, ZETA ) ; % COMPUTE NEW VALUES OF VORTICITY BY ADI E := PI « DX ; LOOP < := 1,1, LL1 DO 3EGIN ZETAIK1 := ZETAIK1 » E * E ; ITY AN 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 090400 090500 090600 090700 090800 090900 091000 091100 091200 091300 091400 091500 091600 091700 091600 091900 092000 092100 092200 092300 092400 092500 092600 092700 092800 092900 093000 093100 093200 093300 093400 093500 093600 093700 093800 093900 094000 094100 094200 094300 094400 094500 094600 094700 094800 094900 095000 095100 095200 095300 095400 095500 095600 095700 095800 095900 096000 096100 096200 096300 81 E := E • DX ; END ; POISSONDIRECT(PSI»ZETA»A»SS»NTSI»NETA»HXtHY»IS»PP) ; E := PI » DX ; LOOP K := 1»1»LL1 DO BEGIN ZETACK] := ZETACK] / ( E * E ) ; E := E * DX ; END ; % COMPUTE NEW VALUES OF THE STRE % BY THE MODIFIED SUCCESSIVE % METHOD IF ST NEQ PINC THEN SIMWRITE (LINE* "NUMBER OF ITERATIONS = »♦ SURFACEVORTICITY ( PSI» ZETA ) ; % ESTIMATES SURFACE VORTICITY BY % THE STREAM FUNCTION VALUES OUTPUTRESULTS (ST»T) ; IF T LSS TLMT THEN GO TO REPEAT ; AM FUNCTION OVERRELAXATION NITER ) ; REFERENCE TO NEAR THE SURFACE END ; % MAIN ; FINi: END . 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 096400 096500 096600 096700 096800 096900 097000 097100 097200 097300 097400 097500 097600 097700 097800 097900 098000 098100 098200 098300 098400 098500 098600 098700 098800 82 xo := NXENO ; YD := nyend ; ZD := NYEND ; IF XD GTR YL) THEN ZD := NXEND BEGIN DOUBLE ARRAY ZETA[0:XDt APPENDIX C % CICHY/MASTER 00000100 BEGIN 00000200 REAL DPSIt TIME2t GOPt YINF, TIMElt REYNt DTIME t STARTT IME ; 00000300 INTEGER NXEND«NYEND»XD»YD»ZD»NITER»LNP»LN I TERt POTENTIAL* 00 0040 SAVE tLN»Pt RECTANGLE tNXLIMltNXLIM2t 00 00 0500 NXLIM3»NYLIM1.NYLIM2»NYLIM3*H0DS» 00000600 PPZtDSVtPZtPP ; 00000700 FILE DATA(MAXRECSIZE=14tRLOCKSI7E=420tKIND=ltACCESS=0t 00000800 TITLE="CICHY/DATAMASTER. M ) I 00000900 READ (DATA »/»REYN»Y INF ♦LNP»LNPP«LNITER»TIME2»DTIME*DPSI»TIME1» 00001000 GOP t NI TER t POTENT I ALt STARTT I ME tSAVEtMGDSt 00001100 DSV t PPZtPZtPPtNXENDtNYENUt RECTANGLE tNXLlMl, 0000120 NXLIM2tNXLlM3tNYLIMltNYLlM2tNYLIM3 ) ; 00001300 00001400 00001500 00001*00 00001700 00001800 D] t PSI(0:XD«0:YD) ♦ 00001900 EElOtYDlt ESlO'.YDlt REIOtYDJ ♦RS[0:YU1» 00002000 RMOSIOrYD) . SMDT[0:YD1 t F[0:ZDlt 00002100 DC0:ZD1 t G(0:ZD)« A[0:5t0:51. B(0:5t0:l)t 00002?00 RZETA[0:ZD)«SLOPEtO:XDl i 00002300 REAL ARRAY ANGLEIOtXD) ; 00002400 ARRAY TITLEARRAYtO: 17] ; 00002500 DOUBLE DPlt DP4, RXSt RYSt RDS4, RERX2t R£RY2t *FBZ2t 00002600 OELXt DELYt Pit DELXSt DELYSt PISt RERXlt ' 00002700 RERX4, RERYlt RERY<4t SIGMAt DELPI t KXYS ; 00002800 REAL ETAt DELlt FKt OELTIt TYMEt Fit WFBZt WFPt 00002900 PROTlMt IOTIM. ELAPTlMt TOTTIM I 00003000 INTEGER ITlt ONE t NPt K» Jt ITZt PROTIMIN. IOTIMINt . 00003100 It Lt NXENDMt NXENDPt NYENOMt NYENDPt 00003200 ELAPTIMlNt TOTTIMSTARTt PROTlMOUTt 00003300 IOTIMOUTt ELAPTlMOUTt PGNUMt HALFt 00003400 NTESTt EVENt NPRESt NISt TOTT iMENDt INP t Q0003500 SECTIONtNPPtNXFILStNXFILFtNYFILSt 00003600 NYFILF,IRECI 00003700 POINTER FTITLE ; 00003800 ALPHA PPRNTt ZPRNTt EPRNTt APRNTt BPRNTt NPRNTt IPRNTJ 00003900 LABEL H0P999t HOPEND ; 00004000 FILE LlNt(KIND=135tMYUSE=2tMAXRECSIZF=22tBUFFERS=2t 00004100 INTM0DE = 3 ) ; -.00004200 FILE CARD(KlND=9t INTMODE=3tMAXRECSIZE=10 ) ; 00004300 FILE STOREP < K IND=DI SK ,F ILETYPE=7t INTMODE=0 tPROTECT ION=0 t 00004400 SAVEFACTOR=99tAREAS=20tAREASIZE=30t ■ * 00004500 TITLE="CICHY/PSI.") ; * 00004600 FILE STOREZ ( K IND=DISK tF lLETYPE=7t INTMGDE=0 tPROTECT ION=0 t 00004700 SAVEFACTOR=99tAREAS=20tAREASIZE=30t 00004800 TITLE="CICHY/ZETA.») ; 00004900 FILE 5TOREC< K IND=DISK tF ILETYPE=7t INTMODE=0 tPROTECT ION=0 t 00005000 SAVEFACT0R=99tAREAS=l t ARE ASI ZE=50 t 000 05100 TITLE = "CICHY/CONSTANTS.'«) ; 00005200 FILE SAVEP ( KIND=DIS»»«<»00047300 00047400 00047500 00047600 00047700 90 NXFILF := NXENOM ; 00047800 NYFILS := 2 i 00047900 NYFILF := NYENDM ; 00048000 ENO ; 00048100 % 00048300 PROCEDURE THOMAS ( NF ILS»NFI|_F, E ♦ TZETA ) ; 00048400 % THE THOMAS METHOD FOR IMVERTING 00048500 % TRIDINGONAL MATRICIES 00048600 DOUBLE ARRAY Et»] » TZETAt*] ; 00048700 INTEGER NFILS»NFILF ; 0OO48R00 BEGIN 00048900 DOUBLE ARRAY WR(0:ZD ]t VR10:ZD ]♦ RRtOtZD ] ; 00049000 INTEGER I» J» K »L . K2 ; 00049100 WRCNFILS) := DINFILS] / EINFILS3 ; 00049200 VRINFILS) := GINFILS] / EtNFILS] ; 00049300 FOR K := NFILS+1 STEP 1 UNTIL NFILF DO 00049400 BEGIN 00049500 RRIK] := EIK) - FtKl « WRIK-l] ; 00049600 KRCK] := DIK] / RR[K1 *» 00049700 END ; 00049R00 FOR K := NFILS+1 STEP 1 UNTIL NFILF DO 00049900 BEGIN 00050000 VRCK] := (GIKJ - F[K)»VR[K-1) ) / RR[K] ; 00050100 END ; 00050200 L := ^FILF ; 00050300 TZETA(L) := VRIL1 ; 00050400 K2 := NFILF -2 ; 00050500 FOR K := NFILS-1 STEP 1 UNTIL K2 DO 00050600 BEGIN 00050700 L := NFILF - K ; 00050800 TZETAtLl := VRIL1 - WRIL1 » TZETACL+11 ; 00050900 END ; 00051000 END ; 00051100 % 00051300 PROCEDURE VORTE ; 00051400 % FIRST HALF OF ADI 00051500 BEGIN 00051600 DOUBLE ARRAY TEIOtXD ] ♦ TZETAtO.'XDl ; 00051700 DOUBLE PFR1» PFR2. Q* S ; 00051800 INTEGER I» J» K »L ; 00051900 FOR I := NXFILS STEP 1 UNTIL NXFILF DO 00052000 BEGIN 00052100 TZETAtll := ZETA(I,NYFILS-11 ; 00052200 END ; 00052300 FOR J := NYFILS STEP 1 UNTIL NYFILF DO 00052400 BEGIN 00052500 FOR I := NXFILS STEP 1 UNTIL NXFILF DO 00052600 BEGIN 00052700 PFR1 := -RDS4 * ( PSIII + liJ) - PSHI-l.J) ) ; 00052800 PFR2 := -RDS4 * ( PSIUfJ+11 - PSIU.J-11 ) ; 00052900 Fill := PFR2 - RERX2 ; 00053000 TEII) := EEtJl ; 00053100 Dill := - PFR2 - RERX2 ; 00053200 Q := RERY2 + PFR1 ; 00053300 S := RERY2 - PFR1 ; 00053400 GUI := Q * RZETAII] + RE(J] ♦ ZETAUtJ] 00053500 + S » ZETA[I»J+1'1 ; 00053600 RZETAII] := ZETAIIfJ] ; 00053700 PROCEDURE VORTS • * % SECOND BEGIN DOUBLE ARRAY TZETAIO: YD) DOUBLE PFRlt PFR2t Q, S INTEGER It Ji » K tL ; 91 ZETA[ItJ-l) := TZETAU1 ; 00053800 END ; 00053900 THOMAS ( NXFILSt NXFILFtTE t TZETA ) ; 00054000 END ; 00054100 FOR i:= NXFILS STEP 1 UNTIL NXFILF DO 00054200 8EGIN 00054300 ZETAIItNYFILFI : = TZETAtll ; 00054400 END ; 00054500 END; 00054600 %«ooo«««»*»«»»««#»*»*»»»»«»««»»»**«»***»*»*«»«w*»»»*«»»#»<»««*#<»<>#o«-tt»»t»00054700 % 00054800 00054900 HALF OF ADI 00055000 00055100 00055200 00055300 00055400 FOR J := NYFILS STEP 1 UNTIL NYFILF DO 00055500 3EGIN 00055600 TZETACJ] := ZETAINXFILS-ltJJ ; 00055700 END; 00055800 FOR I := NXFILS STEP 1 UNTIL NXFILF DO 00055900 BEGIN 00056000 FOR J := NYFILS STEP 1 UNTIL NYFILF DO 00056100 BEGIN 00056200 PFR1 := -RDS4 * (PSHI + ltJl - PSI(I-ltJ) ) ; 00056300 PFR2 := -RDS4 * (PSHItJ+11 - PSICltJ-1) ) ; 00056400 FIJI := -PFR1 - RERY2 ; 00056500 DIJ1 := PFR1 - RERY2 ; 00056600 := RERX2 - PFR2 ; 00056700 S := RERX2 + PFR2 ; 00056800 GUI := Q * RZETAU1 + RSIJ] * ZETAIItJl 00056900 + S * ZETAH + ltJl ; 00057000 RZETACJ] := ZETAUtJl ; 00057100 ZETACI-ltJ] := TZETAIJ] ; 00057200 END ; 00057300 GINXFILS] := GlNXFILS) - F(NXFILS) » ZETAt I t NXF ILS-1 ] ; 00057400 THOMAS ( NYFILSt NYFILFt ES t TZETA ) ; 00057500 END ; 00057600 FOR J := NYFILS STEP 1 UNTIL NYFILF DO 00057700 BEGIN 00057800 ZETAINXFILFtJ) := TZETA[J1 ; 00057900 END ; 00058000 END ; 00058100 PROCEDURE ADI ; 00058300 BEGIN 00058400 LABEL L2 ; 00058500 IF RECTANGLE EQL 1 THEN 00058600 BEGIN 00058700 NETG ; 00058800 FOR I := NXFILS STEP 1 UNTIL NXFILF DO 00058900 BEGIN 00059000 RZETAU] := ZETAt I tNYFILS-1 ) ; 00059100 END ; 00059200 VORTE ; 00059300 FOR J := NYFILS STEP 1 UNTIL NYFILF DO 00059400 BEGIN 00059500 RZETAUJ := ZETAINXFlLS-ltJ) ; 00059600 END ; 00059700 92 vorts ; 00059800 50 TO L2 ; 00059900 end ; 00060000 NETA ; 00060100 FOR I := NXFILS STEP 1 UNTIL NXFILF DO 00060200 3EGIN 00060300 RZETA(I) := ZETAf I ♦NYFIlS-l] ; 00060400 END ; 00060500 vorte ; 00060600 NETB ; 00060700 vorte ; 00060800 netc ; 00060900 vorte; 00061000 NETD ; 00061100 FOR J := NYFILS STEP 1 UNTIL NYFILF DO 00061200 3EGIN 00061300 RZETAIJ] := ZETAfNXFILS-l»J] ; 00061400 END ; 00061500 vorts ; 00061600 NETE ? 00061700 vorts ; 00061800 NETF ; 00061900 vorts ; 00062000 L2: END ; 00062100 ^tttttttttttttt^ttttttttfcttttttattttHs-fttttttttttttttttttttttt^w^tttttttttt^tttttt^tttttttttttttttt % 00062300 PROCEDURE MODSTREAM ; 00062400 BEGIN 00062500 DEFINE NEWPSI = 00062600 3EGIN 00062700 TPSI :=DP1 * PSIIhJ] + DP4 * ( RXS * (PSIII- L»J1 00062800 + PSIt I + 1»J]) + RYS * (PSIt I ♦J-ll + PSI(I»J+1 ]) 00062900 - RMOSt J] * ZETAf It J] ) • 00063000 TP := DABS( TPSI - PSIfltJ] ) ; 00063100 IF TP GTR DPSI THEN GOP := 2.0 ; 00063200 =>SI(ItJl := TPSI ; 00063300 END * X 00063400 INTEGER Ct PASS* SW. I, j; 00063500 DOUBLE TPSI ; 00063600 REAL TP ; 00063700 LABEL H0P444t HOP405t HOP490t Lit L2 ; 00063800 NIS := 1 ; 00063900 HOP405:GOP:=1.0 ; 00064000 SW := ; 00064100 PASS := 1 ; 00064200 H0P444tC := IMP ; 00064300 Li: CASE C OF BEGIN NETD;NETE ;NETF ; GO TO L2 INETG L2;end ; 00064400 C := C + 1 ; 00064500 FOR I := NXFILS STEP 1 UNTIL NXFILF DO 00064600 BEGIN 00064700 FOR J := SW ♦ NYFILS STEP 2 UNTIL NYFILF DO NErfPSI ; 00064800 IF SW = THEN SW:=1 ELSE SW:=0 ; 00064900 END ; 00065000 GO TO LI 1 00065100 L2: IF PASS=1 THEN BEGIN PASS:=2; SW:=i; GO TO H0P444; end ; 00065200 IF G0° LSS 1.5 THEN GO TO H0P490 ; 00065300 NIS := NIS + 1 ; 00065400 ITZ := TIME(2) ; % ELAPSED TIME OF JOB IN 1/60 SECS 00065500 TYME := ITZ / 60 ;%TIME IN SECONDS 00065600 IF TYME LSS TIME2 THEN GO TO HOP405 ; 00065700 93 HOP490:END ; 00065800 %««««»«»««« u »*»»»*#«*'»* »»•»»«*»»•»»■»■»•»»#•»»•»•» #«*«##»»»« it #•»■»•»# #*##»##»*#### o 6 5 9 % 00066000 PROCEDURE SWEEP ; 00066100 BEGIN 00066200 DEFINE NEWPSI = 00066300 BEGIN 00066400 TPSI :=DP1 • PSIUvJ] + DP4 * ( RXS • (PSKI-ltJ] 00066500 ♦ PSItl + l.Jl) + RYS » (PSHI.J-ll + PSI t 1 * J+l ] ) 00066600 - RMOSIJ] * ZETAtlfJ] ) ; 00066700 TP := DABS( TPSI - PSIU.JJ ) ; 00066800 IF TP GTR DPSI THEN GOP := 2.0 ; 00066900 PSI(ItJ) := TPSI ; 00067000 END # ; 00067100 INTEGER I» J; 00067200 DOUBLE TPSI ; 00067300 REAL TP ; 00067400 FOR I := NXFILS STEP 1 UNTIL NXFILF DO 00067500 BEGIN 00067600 TOR J := NYFILS STEP 1 UNTIL NYFILF DO NEWPSI ; 00067700 END; 00067800 END ; 00067900 %»««»«*#«*»»»»»«»#«*«*»»##»*«#»*»»#*«»*»«w<»»«**«»ttw#tttt»»«»*'tt»»»»»»*««»»0O068 % 00068100 PROCEDURE STREAM ; 00068200 BEGIN 00068300 LABEL HOP405* HOP490, Lit L2 ; 00068400 INTEGER C ; 00068500 IF MODS = 1 THEN BEGIN MODSTREAM ; GO TO HOP490 ; END ; 00068600 NIS := 1 ; 00068700 HOP405:GOP := 1.0 ; 00068800 C := INP ; 00068900 Li: CASE C OF BEGIN NETD; NETE ;NETF ; GO TO L2;NETG;G0 TO L2;END ; 00069000 SWEEP ; 00069100 C := C + 1 ; 00069200 GO TO LI ; 00069300 L?: IF GOP LSS 1.5 THEN GO TO HOP490 ; 00069400 NIS := NIS + 1 ; 00069500 ITZ := TIME(2) ; % ELAPSED TIME OF JOB IN 1/60 SECS 00069600 TYME := ITZ / 60 ; % TIME IN SECONDS 00069700 IF TY^E LSS TIME2 THEN GO TO HOP405 ; 00069800 HOP490:END ; 00069900 %«o«»»»««»o*«*»«»«»o««tt««*#*»»»»#««»#tt»»»tt«»»»««»«tt«»«^«*tt-»»»tt«u«-#««<>«-»00070000 % 00070100 PROCEDURE SURFACEVORTICITY ; 00070200 % COMPUTE THE SURFACE VORTICITY 00070300 BEGIN 00070400 INTEGER I* J» K *L ; 00070500 FOR I := 2 STEP 1 UNTIL NXENDM DO 00070600 BEGIN 00070700 ZETAI I»l i: = (8.0 « PSI(I»2J - PSIUOl ) • WFBZ2 ; 00070800 END ; 00070900 END ; 00071000 % 00071200 PROCEDURE SHORTTIMESOLUTION ( T ) ; 00071300 DOUBLE T ; 00071400 BEGIN 00071500 DOUBLE ARRAY StOtXDl ♦ CtOtXDJ* FDZT[0:XD) ; 00071600 DOUBLE VAR1»VAR2»VAR3*VAR4,VAR5,VAR6*ALF»CST1»CST2*CST3»CST4, 00071700 9h Li: INTEG LABEL CST1 CST2 CST3 CST4 CST5 CST6 CST7 CST8 CST9 CST10 CST11 CST12 CST13 CST14 CST15 CST16 CST17 CST18 CST19 CST20 CST21 CST22 CST23 CST24 CST25 CST26 CST27 CST28 CST29 CST30 CST31 ETA : FOR I CST5*C CST13, CST20, INT , VAR8»V VARIO, VAR14, VAR19, ER CC ; L1»L2 ; DSQRT(PI) ; l.o/csTi ; 1.0 / pi ; pi * pi ; CST3 • CST3 ; DSQRT(2.0) ; 1.0 / CST6 ; l.o / 3.0 ; 4.0 « CST2 ; 0.25 * CST1 1.5 » CST1 ; CST8 * CST3 0.5 - 2.0 » 2.0 * CST8 » CST14 * CST1 1.0 + CST15 CST2 * CST16 11.0 » CST2 CST2 » ( CST 8.0 « CST8 * CST2 * ( 8.0 2.0 * T / RE DSQRT( CST22 1.0 / CST23 0.5 « CST24 8.0 » CST8 » REYN /( 8.0 4.0 * CST22 DSQRT(CST28) 4.0 * CST12 2.0 » CST12 = o.o ; STEP 1 UNT ST6,CST7»CST8»CST9,CST10»CST11,CST12» CST14,CST15,CST16»CST17»CST18»CST19, CST21»CST22,CST23»CST24,VAR7»CST25»CST26, TSI, ETA » DTSI» LEVEL* AR9»VAR12»VAR13» CST27* CST28» CST29* CST30, VaRII* PITSI»VAR15»VAR16»VAR17,\/AR18» VAR20»VAR21»VAR22»VAR23,CST31»VAR24 ; ; CST12 ; CST2 ; / 6.0 ; 15 - 1.5 ) ; CST7 * CST2 ; » C5T8 yn ; ) ; ; i CST7 ; » t ) ; * ? * * CST7 - CST15 - 1.0 ) := 1 3EGIN SID cm := DSINfPI* := DCOS(PI* SfNXENDP-I ] := S CCNXENDP-I ] := - ETA := ETA + DEL END ; IL NXENDM/2 DO ETA) ; ETA) ; [11 ; cm ; x ; SCNXE CtNXE TSI DTSI CC : CASE CC : = IF NY FOR J NDP/2J := 1.0 ; MDP/2) := 0.0 ; = 0.0 ? := DELY ; INP 5 cc of begin imeta;metb;netc;go to l2;netg;go to l2;end; cc + i ; F"ILS = 2 THEN NY! := NYFILS STEP BEGIN fils := 1 ; 1 UNTIL NYFILF DO 00071800 00071900 00072000 00072100 00072200 00072300 00072400 00072500 00072600 00072700 00072800 00072900 00073000 00073100 00073200 00073300 00073400 00073500 00073600 00073700 00073800 00073900 00074000 00074100 00074200 00074300 00074400 00074500 00074600 00074700 00074800 00074900 00075000 00075100 00075200 00075300 00075400 00075500 00075600 00075700 00075800 00075900 00076000 00076100 00076200 00076300 00076400 00076500 00076600 00076700 00076800 0007690C 00077000 0007710C 00077200 00077300 00077400 00077500 00077600 00077700 95 ) =>ITSI := PI • TSI ; VARI := DEXP ( PITSI VAR2 := 1.0 / vari ; ALF := CST25 » (VARI - 1.0) ; IF TSI LSS 1.0(5>(«-15 THEN BEGIN ALF := CST25 « ( PITSI * ( 1.0 + 0.5 • PITSI + CST8 * PITSI ) ) ) ; END ; VAR3 := DEXP(-ALF * ALF ) ; VAR4 := ERF( ALF ) ; VAR5 := 1.0 - VAR4 ; VAR6 := ERF( CST6 » ALF ) ; VAR7 := 1.0 - VAR6 ; VAR11 := VAR2 « VAR2 ; VAR12 := CST29 » VAR2 ; VAR13 := 4.0 • CST28 * VAR11 ; VAR14 := VARI - VAR2 ; VAR17 := -VAR2 + VAR3 ; VAR18 := VAR3 * ( 11.0 * VAR5 - - 6.0 ; = VAR3 - 1.0 ; = 1.0 + VAR11 - 2.0 « VAR5 = VAR11 - VAR3 ; = VAR5 - VAR3 ; LSS 1.0i5>K»-15 THEN » ( 1.0 9.0 ) ♦ 4.0 * VAR5 VAR19 VAR21 VAR22 VAR23 IF TSI BEGIN VAR14 VAR15 VAR16 VAR20 VAR24 2.0 ♦ CST8 » PITSI * PITSI ) VAR21 IF ALF BEGIN VAR17 := PITSI • ( := ALF * ALF ; := VAR1S * ( - 1.0 ♦ 0.5 « VAR15 * (1.0 - CST8 * VAR15 ) ) ; := 2.0 * PITSI » := VAR20 * ( -1.0 ♦ 0.5 » VAR20 ♦ ( 1.0 - CST8 * CST20 ) ) ; := VAR24 + 2.0 » VAR4 ; LSS 1.0iap-7 THEN VAR22 END ; := PITSI * (1.0 + CST8 := VAR21 - ( 1.0 - » PITSI VAR16 ; 0. )) ► PITSI VAR16 ; IF END ; ALF LSS BEGIN VAR15 VAR16 1.OP0-7 THEN := ALF • ALF I := VAR15 * ( - 1.0 (1.0 - CST8 * VAR15 + 0.5 > ) ; « VAR15 * VAR18 := 2.0 * VAR16 - VAR4 VAR16 ) ; VAR19 := VAR16 ; VAR23 := - VAR4 - VAR16 ; END ; ( 15.0 - 11.0 » LEVEL := INT VAR14 ♦ CST23 • CST9 * ( VAR17 - CST1 ♦ ALF * VAR5 + CST23 » ( -CST10 * VAR4 + CST11 » ALF * ALF * VAR5 - 1.5 * ALF * VAR3 := ALF * ( CST12 * VAR3 * VAR3 + VAR5 * ( 0.5 * VAR4 - CST31) + ALF » ( VAR3 * (-CST14 * VAR5 CST17) + ALF * VAR5 * (CST8 * VAR5 - CST16) ) ) ) 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 ;ooo 000 ♦ 000 ) 000 77800 77900 78000 78100 78200 78300 78400 78500 78600 78700 78800 78900 79000 79100 79200 79300 79400 79500 79600 79700 79800 79900 80000 80100 80200 80300 80400 80500 80600 80700 80800 80900 81000 81100 81200 81300 81400 81500 81600 81700 81800 81900 82000 82100 82200 82300 82400 82500 82600 82700 82800 82900 83000 83100 83200 83300 83400 83500 83600 83700 96 +CST2 « ( CST26 • VAR6 + CST15 * 00083800 VAR19 + 0.5 * CST8 « VAR18 ) ; 00083900 VAR8 := 8.0 » CST23 * T • INT ; 00084000 VAR9 := VAR5 * < 3.0 + VAR2 * ( 6.0 » CST23 • ALF 00084100 -2.0 ♦ VAR2 * CST23 * ( ALF * ( 4.0 - 00084200 6.0 » CST23 * ALF ) - CST23 ))) + VAR3 * 00084300 2.0 * CST2 » ( -ALF + CST24 - 2.0 * 00084400 CST23 » VAR2 * < 1.0 ♦ VAR2 * ( 1.0 - 00084500 1.5 * CST23 * ALF ))) + CST22 * VAR2 » 00084600 VAR2 ; 00084700 VAR10 := 4.0 » T « CST25 ♦ ( VAR5 » ( VAR12 * 00084800 CST13 - VAR13 * CST14 + ALF * ( -6.0 * 00084900 CST16 - VAR13 » CST13 - VAR3 » VAR12 * 00085000 3.0 * CST2 + ALF * ( -3.0 * VAR12 * 00085100 CST16 + VAR3 * ( 2.0 * CST2 + VAR13 * 00085200 CST14 ) + ALF • ( VAR13 * CST16 - 00085300 VAR5 • VAR13 » CST8 ))) + VAR5 « 00085400 ( -0.5 * VAR12 +ALF » ( 2.0 + 0.5 * 00085500 VAR13 + ALF « VAR12 )) - VAR3 • ( CST2 + 00085600 CST18 * VAR13 )) + VAR3 * ( 2.0 * CST2 « 00085700 ( 1.5 + CST30 ) - CST30 * VAR12 - 00085800 VAR13 » CST19 + VAR3 » VAR12 « 2.0 » 00085900 CST3 + ALF * ( 8.0 » CST12 + CST9 » 00086000 VAR12 * < 1.0 + CST12 ) - VAR3 * ( 2.0 ♦ 00086100 CST3 + CST12 * VAR13 ) - ALF » CST2 * 00086200 ( 2.0 + VAR13 • CST16 ))) + 00086300 VAR13 <* CST20 ♦ VAR7 - VAR13 • CST21 ) ; 00086400 FOR I := NXFILS STEP 1 UNTIL NXFILF 00 00086500 BEGIN 00086600 PSIIIiJ) := ( LEVEL + VAR8 * CC I ] ) « Sill ; 00086700 ZETACI^J] := ( VAR9 + VAR10 » CU) ) * Sill; 00086800 FDZTII1 := PI * ( 8.0 * CST2 * CST29 + 7.0 - 00086900 CST9 * CST25 ♦ REYN » ( CST2 « C5T29 * 00087000 ( 1.0 + CST30 ) - 2.0 ) » Ctl] ) » 00087100 sm ; 00087200 SLOPE! I] := -FDZTCI1 * CST3 ; 00087300 END ; 00087400 TSI := TSI + DTSI ; 00087500 EMD ; 00087600 GO TO LI ; 00087700 L2: TIME1 := T ; 00087800 END ; 00087900 % 00088100 PROCEDURE INITIALIZENET ; 00088200 % INITIALIZE THE PSI NET WITH THE POTENTIAL FLOW00088300 % SOLUTION* RELAX IT* COMPUTE THE SURFACE 00088400 % VORTICITY»AND SET THE REST OF THE VORT ICI TY00086500 % NET EQUAL TO ZERO 00088600 BEGIN 00088700 INTEGER NXP» NXM, II * I»J,K»L : 00088800 REAL TCG» ETCG, RETCGS, FTG, DIE* SI 1 • PIYIN. CSIGM, 00088900 PIDEL* PIE ; 00089000 LABEL LI ; 00089100 IF NITER NEQ THEN BEGIN FETCHPSIZETA ; GO TO LlfENDl 00089200 NXM := NXEND / 2 ; 00089300 NXP := NXM ♦ 1 ; 00089400 SIGMA := 0.0 ; 00089500 FOR J := 1 STEP 1 UNTIL NYENDM DO 00089600 BEGIN 00089700 97 ' TCG := PI « SIGMA ; ETCG := EXP(TCG) ; RETCGS := 1.0 / (ETCG * ETCG) ; FTG := ETCG * (1.0 - RETCGS) ; ETA := 0.0 ; FOR I := 1 STEP 1 UNTIL NXP DO BEGIN PIE := PI • ETA ; SI1 := SIN(PIE) ; psn itji := ftg » sii ; ETA := ETA ♦ DELX ; END ; sigma := sigma ♦ dely ; end; piyin := pi * yinf ; csigm := exp(piyin) ; ftg := csigm ; if potentials then BEGIN RETCGS := 1 .0/ (CSIGM»CSIGM) ; FTG := CSIGM ♦ (1.0 - RETCGS) I END ; PIDEL := PI * DELX ; FOR I := 1 STEP 1 UNTIL NXP DO BEGIN Fl := i-i ; »si[i» nyend1 := ftg » sin(pidel » fi) ; end; for j := 1 step 1 until nyend do BEGIN FOR I := 1 STEP 1 UNTIL NXM DO BEGIN II ;= NXENDP - I ; PSIIII. Jl := PSltl.Jl ; end ; end ; for j := 1 step 1 until nyend do BEGIN p"OR I := 1 STEP 1 UNTIL NXEND DO BEGIN ZETAU»J] := 0.0 ; END ; END ; IF PP=1 THEN PRINT1 ( PSI » IPRNT ) ; IF STARTTIME EQL 0.0 THEN BEGIN STREAM; SURFACEVORT LI* IF RECTANGLE NEQ 1 THEN INP := ; IF STARTTIME GTR 1.0ia«-30 AND SECTION = 1 THEN SHORTT IMESOLUT ION ( STARTT IME) END ; % PROCEDURE MATINV ( A»N»B»M, DETERM ) ; DOUBLE ARRAY At»»*l ♦ B[*»») ; INTEGER N»M ; DOUBLE DETERM ; BEGIN DEFINE IROW = JROW #♦ ICOLUM SWAP = T # ; DOUBLE ARRAY PIVOT(0:5] ; DOUBLE T ; = JCOLUM n% AMAX = T #♦ 00089800 00089900 00090000 00090100 00090200 00090300 00090400 00090500 00090600 00090700 00090800 00090900 00091000 00091100 00091200 00091300 00091400 00091500 00091600 00091700 00091800 00091900 00092000 00092100 00092200 00092300 00092400 0009250C 00092600 00092700 00092800 00092900 00093000 00093100 00093200 00093300 00093400 00093500 00093600 00093700 00093800 00093900 00094000 00094100 icity; end ;ooo9420o 00094300 00094400 ; 00094500 00094600 «««»<»«««««#* 00094700 00094800 00094900 00095000 00095100 00095200 00095300 00095400 00095500 00095600 00095700 98 INTEGER ARRAY IPI VOTl : 5 ) * INDEXI : 5, : 2 ] ; 00095800 INTEGER JROW, JCOLUM, LI »I»J»K»L ? 00095900 LABEL HOP720 ♦ HOP777 ; 00096000 DETERS := 1.0 ; 00096100 FOR J := 1 STEP 1 UNTIL N DO 00096200 BEGIN 00096300 IPIVOTtJ) := : 00096400 END ; 00096S00 FOR I := 1 STEP 1 UNTIL N DO 00096600 BEGIN 00096700 LABEL HOP260, HOP380< HOP550 ; 00096800 AMAX := 0.0 ; 00096900 FOR J := 1 STEP 1 UNTIL N DO 00097000 BEGIN 00097100 LABEL HOP105 ; 00097200 IF (IPIVOTl Jl-1) EQL THEN GO TO HOP105 ; 00097300 FOR K := 1 STEP 1 UNTIL N DO 00097400 BEGIN 00097500 LABEL HOP100 ; 00097600 IF (IPIVOTtK) -1) GTR THEN GO TO HOP777 ; 00097700 IF (IPIVOTtK) -1) EQL THEN GO TO HOP100 ; 00097800 IF (DABS(AMAX) - DABS(AIJ,K) )) GTR 0.0 THEN 00097900 GO TO HOP100 ; 00098000 00098100 00098200 00098300 HOP100: END ? 00098400 HOP105: END ; 00098500 00098600 00098700 00098800 FOR L := 1 STEP 1 UNTIL N DO 00098900 BEGIN 00099000 SWAP := AlIROWfL) ? 00099100 AtIROW,L) := AtICOLUM,L) ; 00099200 At ICOLUM, L) := SWAP ; 00099300 END ; 00099400 IF M LEQ THEN GO TO HOP260 ; 00099500 r OR L := 1 STEP 1 UNTIL m DO 00099600 BEGIN 00099700 SWAP := BlIROW,L 1 ; 00099800 BtIROW,L) := BtICOLUM,L) 5 00099900 BtICOLUM,L] := SWAP ; 00100000 end; ooiooioo hop260: indexu,1) := irow ; 00100200 INDEXtI,2) := ICOLUM ; 00100300 ?IVOTtI) := At ICOLUM, ICOLUM] ; 00100400 DETERM := DETERM » PIVOTtI] ; 00100500 IF PIVOTtI) EQL THEN GO TO HOP720 : 00100600 At ICOLUM, ICOLUM) := 1.0 ; 00100700 FOR L := 1 STEP 1 UNTIL N DO 00100800 BEGIN 00100900 AtICOLUM,L) := AtICOLUM,L) / PIVOTtI) ; 00101000 END ; 00101100 IF M LEO THEN GO TO H0P380 ; 00101200 FOR L := 1 STEP 1 UNTIL M DO 00101300 BEGIN 00101400 BtICOLUM,L) := BtICOLUM,L) / PIVOTtI) ; 00101500 END ; 00101600 HOP380: FOR LI := 1 STEP 1 UNTIL N DO 00101700 IROW := j ; ICOLUM : = k ; AMAX := AtJ»K] ; END ; END • IP IVOTtlCOLUM] := IPIVOTl ICOLUM) + l ; IF DE (IROW TERM := - ICOLUM) : - DETERS EQL ; THEN GO TO HOP260 99 BEGIN IF (Ll-ICOLUM) EQL THEN GO TO HOP550 ; T := A[Ll»ICOLUM) ; AILltlCOLUM] := 0.0 ; FOR L := 1 STEP 1 UNTIL N DO BEGIN AtLl.Ll := ACL1»L1 - AUCOLUM.Ll » T ; END ; IF M LEQ THEN GO TO HOP550 ; FOR L := 1 STEP 1 UNTIL M DO BEGIN BIL1.LJ := BtLl.L] - BUCOLUM.L1 »T ; END ; end; HOP550: END ; FOR I := 1 STEP 1 UNTIL N DO BEGIN LABEL HOP710 ; L := N+ 1 - I ; IF 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 101800 101900 102000 102100 102200 102300 102400 102500 102600 102700 102800 102900 103000 103100 103200 103300 103400 103500 103600 103700 103800 103900 104000 104100 104200 104300 104400 104500 104600 104700 104800 104900 105000 105100 105200 105300 105400 105500 105600 105700 105800 105900 106000 106100 106200 106300 106400 106500 106600 106700 106800 106900 107000 107100 107200 107300 107400 107500 107600 107700 100 9EGIN 00107800 FOR J := 1 STEP 1 UNTIL M DO 00107900 BEGIN 00108000 All, J] := 0.0 ; 00108100 FOR K := 1 STEP 1 UNTIL N DO 00108200 BEGIN 00108300 A[I»J] := A[I»J] + CtK»Il»ClKtJ) » WlK] ; 00108400 EN ; 00108500 END ; 00108600 r ND ; 00108700 FOR I := 1 STEP 1 UNTIL M DO 00108800 9EGIN 00108900 311,11 := 0.0 ; 00109000 FOR K := 1 STEP 1 UNTIL N DO 00109100 BEGIN 00109200 BCI*1J := BUtl) + C(K»I) » YlKtl)»W(K] ; 00109300 END ; 00109400 END • 00109500 MATINV ( A,M,B*L, DETERM ) ; 00109600 END ; 00109700 *«. 9«««««a««o««» ^ <>tt ^ w< ). i »tt*#«*«-*tt*«-«<>«^» * ( K IND=D I SK*FILETYPE = 7*INTM0DES* PROTECTIONS* FILE FILE FILE STOREZ ( FILE STOREC( FILE SAVEP 00000100 00000200 00000300 00000400 00000500 oooocsoo 00000700 00000800 00000900 00001000 00001100 00001200 00001300 00001SC 00001500 00001800 00001700 00001800 00001900 00002000 00002100 00002200 00002300 00002400 00002500 00002600 00002700 00002800 00002900 00003000 00003100 00003200 00003300 00003400 00003500 00003600 00003700 00003800 00003900 00004000 00004100 00004200 00004300 00004400 00004500 00004600 00004700 00004800 00004900 00005000 00005100 00005200 00005300 00005400 00005500 00005600 00005700 00005800 00005900 00006000 106 FILE SAVEZ ( FORMAT FORMAT FORMAT FORMAT FORMAT FMT1 ( FMT2 ( FMT3 ( FMT4 ( FMT5 ( FORMAT FORMAT FORMAT FORMAT FORMAT FORMAT FMT6 ( FMT7 ( FMT8 ( FMT16 FMT75 FMT90 ACCESS= MAXRECS KIND=DIS ACCESS= MAXRECS F10.0»?I 6F10.0) 2F10.3* X56 X49 X49 X49 X49 X49 X49 X49 X49 X49 X49 X49 X49 X49 ii "PA "REY "LNP "LNI ••TIM ••TIM •'NIT "POT "NXE "IOX "DEL "STT "DSV "IRE •) ; 1*SAVEF IZE=2»Y K»FILET 1»SAVEF IZE=2*Y lOtFlO. • 3110 ) RAMETER N TER E? El ER ENT MD X IME ACTOR D + 2,T YPE = 7 ACTOR D + 2«T 0) ; = '* —ii = n = ■< = •• = »• = n = •• = •• = m = n = i* — ii S"/) F13.7 15 15 F13.7 F 1 3 . 7 15 15 15 15 F13.7 F13.7 15 15 / 99, TLE="CICHY/SPSI.") ; INTMODE=0 . PROTECT I ON=0* 99* TLE="CICHY/SZtTA.") ; * 102 * 103 * 105 *, 107 X49,"YINF ="*K13.7/ X49*"LNPP =". 15 / X49»"DTIME ="*F13.7// / X49»",\JYEND = "» 15 / X49«"IQY = "» 15 / X49»"DELY = '%F13.7// X49«"SAVE = ", IS // X49»"PPZ ="♦ 15 / X49,"SECTI0N=", 15 ) 5 % 108 % PROCEDURE DATIME ; BEGIN VALUE ARRAY A( "TIME ="» F10.2) : 10F8.4) ; ( //// X ( // ) ; « * -» •tnn* tt to a to to to to a # to to it to to to to to ■::• to to tt » a » to to w a to to a an it a to a- to to a to a >7»"CICHY'SX108, " PAGE "♦ 13// ) ? 60" SUN' 60" OWT: MON"» 60" TUES"»60"WEDNES"» 60" THURS", 60" FRI", 60" SATUR") ; MlN»G»PtY ♦ TIMEON ? (l) ; 100 + 10»Y.[ 11:6] + Y.{5:6] ; 10 + Y. [23:6] ; Y MOD 4=0 THEN 29 ELSE 2*«31t30.31«30»31«31»30.31»30 UO H THEN GO TO OWT ELSE BEGIN D : =D-H : m : =m+ 1 END : 216000 ; IV 3600) MOD 60 ; n begin q:=m+h;p:=y-1 END ELSE HEGIN q:=m-1;p:=Y END ? R(TITLFAPPAY,6) ; RY 6" " FOP 132 : +85 *Y POINTEP(A[ ( (Q«26-2) OlV 10+D+P+P.(li:i0]+1) MOU 7]«6) FOR 6» 6"DAY. "♦ M FOR 2 DIGITS»6"/"» D FOP 2 DIGITS* 6"/'S Y FOR 2 DIGITS* 6"» " « ( 1 2#PEAL (O : =H MOD 12=0)+G FOR 2 DIGITS* 6": M »MIN FOP 2 0IGITS.6" " » IF H GEO. 1? THFN 60"P" ELSE 60"A" FOR 1* 6"M" END datime ; 9&ttttttatt«'il-tttt.««'tttttttt"<}ti ; % PRINTOUT ROUTINE DOUBLE ARRAY NETTYPEI***] ; ALPHA NETNAME ; BEGIN ALPHA DESIG ; INTEGER NC» II» Ml* M2 * KK ♦ NXMARK* NXLIM* SECHAF ; LABEL H0P512* HOP503 * HOP504 ; FORMAT FMT10< X20*"NITER ="I5*X23* "SECTION NUMBER X5* A4, X19* "TIME = "* F7.3/ ) ; FORMAT FMT11 ( 11E12.4 ) ; NC := NYEND / 11 + 1 ; FOR II := 1 STEP 1 UNTIL NC DO BEGIN SECHAF := ; Ml := 11 * II - 10 ; IF (Ml - NYEND) GTR 0.0 THEN GO TO H0P512 ; M2 := Ml + 10 ; IF (M2 - NYEND) LEO 0.0 THEN GO TO HOP504 ; M2 := NYEND ; HOP504: IF HALF=1 THEN BEGIN NXMARK : =NXFND/2 . ; NXLIM:= DESIG := APRNT: END ELSE BEGIN NXMARK : vNXENDl NIXLIM:= 1 5 DESIG := NPRNT ; FND : HOP503: WRITE (LINE [SKIP 11); PGNUM := PGNUM + 1 ; WRITE (LINE* FMT75* PGNUM ) ; WRITE ( LINE* 16 » TITLEARRAYf*) ) : WRITE (LINE* FMT90 ) ; WRITE (LINE*FMT10* NITER* IN DEbIG* NETNAME* 11 FOR I := NXLIM STEP 1 UNTIL NXMARK DO BEGIN WRITE(LINE* FMT11* FOR J := Ml STEP 1 UNTIL NETTYPEt I*J]) ; END ; IF SECHAF = 1 THEN GO TO H0P512 ; IF HALF=1 THEN BEGIN SECHAF:=l; NXLIM:= NXMARK + NXMARK:=NXEND; DESIG := 3PRNT : GO TO HOP503* END ; H0P512: END ; END ; %■»■»»■»»•»»•»« «e*tn>o* a (nt«in»>n>tt# it mm «## # innnnnnnn*inma# * i> it it it a it Hit* % PROCEDURE TYMIN: * INITIALIZE THE TIMING ROUTINE BEGIN PR0TIMIN:= TIME(2) » IOTIMIN := TIMEO) ; ELAPTIMIN := TIME(l) ; END ; % it a * «• a •» it •» a it a it it « «• a a it it a- it it «• a it « it <- it w it » a v.- it it it it it it a a mt n it it a- it it it it « •> a it it it a "*000 000 000 000 000 12100 12?00 12300 12400 12500 12600 12700 12800 12900 13000 13100 13200 13300 13A00 13500 13600 13700 13800 13900 14000 14100 14200 14300 14400 14500 14600 14700 14*00 14900 15000 15100 15200 15300 15<+00 15500 15600 15700 15800 15900 16000 16100 16200 16300 16400 16500 16600 16700 16800 16900 17000 17100 17200 17300 17400 17500 17600 17700 17800 17900 18000 108 PR0TIM0UT := TIME<2) 5 00018100 IOTIMOUT := TIME(3) ? 00018200 ELAPTIMOUT := TIME(l) : 00018300 PROTI^ := (PROTIMOUT - PROT I M I n ) /t>0 ; 00018400 IOTIM := (IOTIMOUT - IOTIMIN)/60 ; 00018SOO ELAPTIM := (ELAPTIMOUT - ELAPTIMIN) /60 ; 00018600 TOTTIMEND := TIME(l) ; 0001*700 TOTTI^ := (TOTTIMEND - TOTT IMST ART ) /3foU0 5 00018*00 WRITE(LINE»FMT51) ; 00018900 WRITE* LINE* FMT50* PROTIM, lOTIMt ELAPTIM, TOTTIM); 00019000 END ; 00019100 %#«#*■»#■»•»*•*«•»*#» tt»tt##tttttttt«tttttttttt#tt#tt»iM*ttttW#tttttt #####««•«•«■',>•«■•& -aw 00019200 % 00019300 PROCEDURE CHECKTIMELIMIT ; 00019^00 BEGIN 00019500 % ELAPSED TIME OF JOB IN 1/60 SECS 00019600 ITZ := TIME<2> ; 00019700 % ELAPSED TIME IN SECS 00019800 TYME := ITZ / 60 ; 00019900 % COMPARE TO TIME LIMIT 00020000 IF TYME GTR TIME2 THEN GO TO H0?999 : 00020100 END J 00020200 % 00020400 PROCEDURE SAVEPSIZETA ; 00020500 BEGIN 00020600 INTEGER RECORD ; 00020700 IF IREC=1 THEN 0002080C BEGIN 00020900 SAVEP.AREASIZE := XO ; 00021000 SAVEP. AREAS := 100 ; 00021100 00021200 WRITE (SAVEPtOl t»«REYN»Y INF* 00021300 DTIME»DPSI»P0IENTIAL« 00021^00 WFP,NXEND,NYEND*DELX»0ELY,STARTT1ME» 00021500 M00S*RECTANGLE*NXLIM1,NXLIM2«NXLIM3» 00021600 NYLIM1»NYLIM2»NYLIM3»TITLEARRAY[«J) 5 00021700 SAVEZ.AREASIZE := XO 5 00021800 SAVEZ. AREAS := 100 : 00021900 WRITE (SAVEZtO]»*«REYN»YINF» 00022000 OTIME^D^SI. POTENTIAL* 00022100 WFP,NXENU*NYEND»DELX»DELY»STARTTIME» 00022200 M00S»RECTANGLE«NXLIM1»NXLIM2«NXLIM3» 0022300 NYLIM1«NYLIM2*NYLIM3»TITL£ARRAY[«J) ; 00022400 LOCK(SAVEP) ; 00022500 LOCK(SAVEZ) ; 00022600 SAVEP.AREASIZE : = ; 00022700 SAVEP. AREAS := ? 00022800 SAVEZ.AREASIZE := ; 00022900 SAVEZ. AREAS := ; 00023000 END ; 00023100 I := 1 ; 00023200 FOR RECORD:=IREC STEP 1 UNTIL XD-1+IREC DO 00023300 3EGIN 00023400 WRITE(SAVEPl RECORD N2*Y0 + 2.PSin»»J) ? 00023500 WR I TE( SAVEZ [ RECORD ]»?*YD + 2tZET A [ I»*J) ! 00023600 I := I + 1 : 00023700 END ; 00023800 IREC := IREC ♦ XO ; 00023900 CLOSE(SAVEP) ; 00024000 109 CLOSE(SAVEZ) ; 00024100 END ; 00024200 % 00024400 PROCEDURE STOREPSIZETA ; 00024500 BEGIN 00024600 STORE . MAXRECSIZE:=2*YD+2 ; 00024700 STORE 3 . AREAS:=20 ; 00024800 STORE . AREASIZE:=JO ; 00024900 FOR l:=l STEP 1 UNTIL XD DO */R I TE (STORED 2»YD+2.PSI I I ♦ » ] ) ; 00025000 LOCK(STOREP) ; 00025100 STOREZ.MAXRECSIZE:=2»YD+2 ; 00025200 ST0REZ.AREAS:=20 ; 00025300 STORE7.AREASIZE':=30 ; 00025400 FOR l:=l STEP 1 UNTIL XD DO rfR I TE ( STOREZ . 2»YD+2. ZETa { I . * ] ) ; 00025500 LOCK(STOREZ) ; 00025600 END : 00025700 % •»■»•»»•»«• tt «<* ■ttttW «• tt « «• ■» «• V tt i> « » » '* -M- «■ tt «• «• ■» * tt « tt it tt ■»■»«•H>H»«^^y»U^tt«y»UttW<>^W1>Wv-H>«^Wtt■&«•«•» J > *- » <* a » * # # « w -K- ^ J > w «■ » * « -;* J ' J > hj- # # # * tt « «■ [EOF21 ; 00031700 E0F2: l:= I-lt 00031800 j:= (STOREP.MAXRECSIZE) UIV 2 X 00031900 CLOSE (STOREP) : 00032000 CLOSE (STOREZ) ? 00032100 ENID ; 0003220C END : 00032300 %tototototototo*tototototototototo Hi- to to •:* « •:<• to to to to a a to •::• •» *- <;• to «• «• w ■«■ W w -a <* tt tt « to <> -i> o <> «• •» -a- •> •> o o 3 2 4 % 00032500 PROCEDURE 'ETCHCONSTANTS t 00032600 BEGIN 00032700 STOREC. AREAS := ; 00032800 STOREC.AREASIZE := 5 0003290C IF STOREC. PRESENT THEN 00033000 3EGIN 00033100 READ(STORFC« * * REYN.YINF,OTI*lE.DPSI *TIMEl»bOP« 00033200 N I TEh 1 . NX END. N YEND ♦ SECT I ON » RECTANGLE* 0033300 NXLIM1,\JXLIM2«NXLIM3»NYLIM1,NYLIM2. 00033400 NYLIM3.M00S»IQX»1QY •IREC) ; 00033500 CLOSE (STOREC) ; 00033600 END ; 00033700 END ; 00033800 fytotototototototototototototoittotototo •:!• to to to « tototo # i> » a «. to to # <:■ to *- to to a to <<• <> -:!• to to to <;• -;:• to *• «• » u * ^ <> *«-><<■ «• «• a «■ <}• v ■> » •» <<• «■ « ■«• « <)■ «■ # ««• ■» <* » ■**•«•«••»>•»«•*•«•* ^00036000 Ill PROCEDURE SETF67 (IQ) ; INTEGER 10 » BEGIN INTEGER N5»I»K1,K»II_; LABEL L11»L10»L9 ; N3 := 2**10 ; N7 := N3/2 ; N5 := N3/4 ; I := 1 ; INOEX[ I] := N5 ; SKI) := 1.0 / SQRK2.0) K := I ; I := 1+1 ; % Lll: IL := I ; IF I EQL N7 THEN GO TO L9 ; L10: Kl := INDEX[K]/2 ; INDEXED := Kl ; Sit I] := DSIN(PI»K1/N3) 5 Kl := N7 - Kl ; I := 1+1 ; INOEXt I ] := Kl ; SIII] := DSIN(PI»K1/N3) ; K := <+l ; I := l+l ; IF K NEO IL THEN GO TO L10 ; GO TO Lll ; L9: END ; %»-»«««»«*»<>«*tt»w^»w^^<>»*«w«««^'>«H>^^**<>»()0052400 00052500 00052600 00052700 00052800 00052900 00053000 00053100 00053200 00053300 RPfK] L + ll 114 TZETAtl] := ZETAU.1J ; 00053400 END ; 00053500 FOR J := 2 STEP 1 UNTIL NYENDM DO 00053600 3EGIN 00053700 ^OR I := 2 STEP 1 UNTIL NXENDM DO 00053800 BEGIN 00053900 PFR1 := -RDS4 » ( PSIII + 1»JJ - PSIU-l.JJ ) ; 00054000 PFR2 := -RDS4 » ( PSItI»J+l] - PSHI.J-1] ) ; 00054100 Fill := PFR2 - RERX2 : 00054200 TEt II := EEt J] : 00054300 Dill := - PFR2 - RERX2 ; 00054400 Q := RERY2 + PFRl ; 00054500 S := RERY2 - PFR1 ; 00054600 G(I) := Q * ZETA(I»J-1] + PECJ] » ZETAU»J) 00054700 + S * ZETA[ItJ+l) ; 00054800 ZETA[I,J-1] := TZETAII] ; 00054900 END ; 00055000 THOMAS ( NXENDM, TE ♦ TZETA ) S 00055100 END; 00055200 FOR i:= 2 STEP 1 UNTIL NXENDM DO 00055300 3EGIN 00055400 ZETA[ I, NYENDM] := TZETAII] ; 00055500 END ; 00055600 END: 00055700 %«««««#»««#«»»«» w ««««««««'»««««««##«# ^««»###»«^««tttttt «««»««« »«tftf#«»««u»tf-»o0055800 % 00055900 PROCEDURF VORTS : 00056000 % SECOND HALF OF ADI 00056100 8EGIN 00056200 DOUBLE ARRAY TZETACOtYDJ ; 00056300 DOUBLE PFRl, PFR2» Q» S ; 00056400 INTEGER I« J» K »L ; 00056500 FOR J := 2 STEP 1 UNTIL NYENDM DO 00056600 BEGIN 00056700 TZETAfJ] := ZETAU»J] ; 00056800 END; 00056900 FOR I := 2 STEP 1 UNTIL NXENDM DO 00057000 3EGIN 00057100 FOR J := 2 STEP 1 UNTIL NYENDM DO 00057200 BEGIN 00057300 PFR1 := -RDS4 « (PSItI+l»J] - PSIII-1»J] ) ; 00057400 PFR2 := -R0S4 * (PSIUfJ+1) - PSIU,J-1) ) ; 00057500 F[J] := -PFRl - RERY2 ; 00057600 DfJ] := PFRl - PEWY2 ; 00057700 := RERX2 - PFR2 ; 00057800 S := RERX2 + PFR2 : 00057900 GUI := U » ZETA(I-1,J] + PStJ] « ZETA£IfJ] 00058000 ♦ S * ZETAtI+l»J] : 00058100 ZETAfI-l,J] := TZETACJ] ; 00058200 END ; 00058300 3121 := GI2] - F[2] » ZETAfI,l] ; 00058400 THOMAS ( NYENDM, ES « TZETA ) : 00058500 END : 00058600 FOR J := 2 STEP 1 UNTIL NYENDM DO 00058700 3EGIN 00058800 ZETAfNXENDM,J] := TZETAfJ] ; 00058900 END ; 00059000 END ; 00059100 %tt««)nnt 0061000 % 00061100 PROCEDURE ZERO(L) ? 00061200 INTEGER L ; 00061300 BEGIN 00061400 INTEGER I ; 00061500 FOR I := 1 STEP 1 UNTIL N2 DO 00061600 BEGIN 00061700 ZtL+I-11 := 0.0 ; 00061800 END ; 00061900 END ; 00062000 %»»«*'U-^-«-«tt 00062100 % 00062200 PROCEDURE ^tt «■»•»«• ^-fjttttw-tt^wttwtt^tttttttt 00072900 % 00073000 PROCEDURE MEG(I1»I3»I2) ; 00073100 INTEGER I1»I3»I2 ; 00073200 BEGIN 00073300 INTEGER K ; 00073400 FOR K:=I1 STEP 12 UNTIL 13 DO 00073500 BEGIN 00073600 YtK] := -YCK] ; 00073700 END ; 00073800 END : 00073900 %»^o»««»»n»#»»4>»««*»««^«^^««nutty«^««^tti*w«^ttWH>tti>wtt^««<»«<>«tt«w«*»««^uwtt« 00074000 % 00074100 PROCEDURE F0UR67(IQ) ; 00074200 INTEGER IQ ; 00074300 BEGIN 0007^400 INTEGER N5»N11»N31»I»IP»JF,I1,I2»I3 : 00074500 N4 := 2»»IQ ; 00074600 N3 := N4 ; 00074700 N5 := N3/4 ; 00074800 N7 := N3/2 ; 00074900 Nil := 3»N7 ; 00075000 N31 := N3+1 : 00075100 ZCN31] := 0.0 ; 00075200 Zt 1) := 0.0 ; 00075300 N2 := N3 ; 00075400 FOR lis 2 STEP 1 UNTIL IQ DO 00075500 3EGIN 00075600 TFOLDUtltZ) ♦ 00075700 ^2 := N2/2 ; 00075800 END ; 00075900 YIN7+11 := ZI21 ; 00076000 JF := N5 I 00076100 FOR I=> := 2 STEP 1 UNTIL 10 DO 00076200 BEGIN 00076300 118 LI := N2+1 ; 00076400 ZERO(l) ; 00076500 ISL := 1 ; 00076600 KFOLD ; 00076700 11 := 3*JF+1 ; 00076800 12 := 4«JF ; 00076900 13 := Il+(N2/2 -1>»I2 ; 00077000 MEG(I1»I3«I2) ; 00077100 M2 := N2+N2 ; 00077200 JF := JF/2 ; 00077300 END ; 00077400 END ; 00077500 %'>ft^»« +F2*K0REl I«J] 00079100 +B*(KOREt I»J+l)+KORE( I • J- 1 )) ; 00079200 K0RE(I-1»J1 := X ; 00079300 X := W ; 00079400 END ; 00079500 «f := 0.0 ; 00079600 W := -Fl» j >««*^<>u^«^tt»*i>««wtttt»^*ui>»<>«^«««^»-u«tt^*i*^w^^w«^««tt«^^^^^-^-»'tJ*^tt-tt-(>'t>-u--tt<»-'Utt 00081400 % 00081500 PROCEDURE STOREX(J) ; 00081600 INTEGER J ; 00081700 BEGIN 00081800 INTEGER I ; 00081900 FOR l:= 2 STEP 1 UNTIL NXENDM DO 00082000 3EGIN 00082100 L25 ; 00083100 DOUBLE ARRAY BBIOtll) ; 00083200 INTEGER N1»NN»J»NV.NW»NU»IP ; 00083300 DOUBLE B ; 00083400 IP := IP1 ; 00083500 BB[ 11 := A ; 00083600 B := A ; 00083700 N4 := ; 00083800 IP := IP-1 ; 00083900 FOR Nl :=1 STEP 1 UNTIL IP DO 00084000 BEGIN 00084100 M4 := N4 + 1 ; 00084200 \|W := 2«»N1 5 00084300 \|V := 2»NW ; 00084400 IF (NV+1) GTR (NYEND-NV) THEN GO TO L21 ; 00084500 FOR J :=(NV+1) STEP NV UNTIL (NYEND-NV) DO 00084600 BEGIN 00084700 KORECItJ] := B*KOREl I»J]+KOREl I ♦ J-NW 1+KORE I I*J+Nw) ; 00084800 END ; 00084900 L2l: 3 := B»B - 2.0 ; 00085000 3BIN1+1] := 8 ; 00085100 IF B LSS 1.068*14 THEN GO TO L2 i 00085200 IF (NV+1) GTR (NYEND-NV) THEN GO TO L12 ; 00085300 FOR J:=(NV+1) STEP NV UNTIL (NYEND-NV) DO 00085400 BEGIN 00085500 KOREU»J) := -KOREC I»J)/B : 00085600 END 5 00085700 GO TO L12 ; 00085800 L2: END ; 00085900 J := (NYENO+D/2 ; 00086000 KOREUfJl := -KORE(I»J]/B ; 00086100 L12: FOR NN:=1 STEP 1 UNTIL N4 DO 00086200 3EGIN 00086300 \U := N4-NN ; 00086400 9 := BBIN1+1 ] ♦ 00086500 NJW := 2**N1 ; 00086600 NJV := 2*NW ; 00086700 MU := 3»NV J 00086800 J := NV+1 ; 00086900 <0REII»J1 := (KOREt I»J+NV] - KORElItJ] ) /B ; 00087000 J := NYEND-NV ; 00087100 *«^*^tt«*tt«««<»««NN»I»NV»NW»NU»IP ; 00088800 DOUBLE B ; 00088900 IP := IP1 ; 00089000 BBtl] := A ; 00089100 00089200 00089300 ; 00089400 STEP 1 UNTIL IP DO 00089500 00089600 N4 + 1 ; 00089700 2»»(N1-1) ; 00039800 2*NW ; 00089900 IF (NV+1) GTR (NXEND-NV) THEN GO TO L21 ; 00090000 TOR I :=(NV+1) STEP NV UNTIL (NXEND-NV) DO 00090100 BEGIN 00090200 KORE[I*JJ := B*KORE[ I*J]+KOREt I-NW ♦ J ] +KORE t I+NW.J] ; 00090300 END ; 00090400 L2l: 3 := B»B - 2.0 ; 00090500 3BIN1+1 J := B ; 00090600 IF B LSS 1.0L«eai4 THEN GO TO L2 ; 00090700 IF (NV+1) GTR (NXEND-NV) THEN GO TO L12 ; 00090800 F"OR l: = (NV+l) STE» NV UNTIL (NXEND-NV) DO 00090900 BEGIN 00091000 KORElItJJ := -K0RE(I»J1/B ; 00091100 END ; 00091200 GO TO L12 ; 00091300 L2: END 5 00091400 I := (NXEND+D/2 ; 00091500 KOREtI*J] := -KORE[I»J]/B ; 00091600 L12: FOR NN:=1 STEP 1 UNTIL N4 DO 00091700 3EGIN 00091800 sjl := N4-NN ; 00091900 3 := B8CN1+1 ] ; 00092000 NV := 2»<*N1 ; 00092100 MU := 3«NV ; 00092200 I := NV+1 ; 00092300 = 2.0 • T / REYN ; CST23 ! = DSQRTC CST22 ) ; CST24 : = 1.0 / CST23 ; CST25 ! = 0.5 » CST24 ; CST26 ! = 8.0 » CST8 * CST7 ; CST27 : = REYN /( 8.0 » T ) ; CST28 i s 4.0 » CST22 ; CST29 ! s DSQRT(CST28) ; CST30 : = 4.0 * CST12 ; CST31 : = 2.0 » CST12 ; ETA : = 0, ,0 ; FOR I : = 1 STEP 1 UNTIL NXENDM/2 DO 9E :gin SI i 1 := DSIN(PI*ETA) ; CI i 1 := DCOS(PI»ETA) ; CST15 - 1.0 ) I StNXENDP-Il := Sill CtNXENDP-Il := -C[I1 ETA := ETA + DELX ; END ; StNXEMDP/21 := 1.0 • CCNXEMDP/2) := 0.0 • TSI := 0.0 • DTSI := DELY ; FOR J := 1 STEP 1 UNTIL NYENDM 3EGIN =>ITSI := PI * tsi ; VARl : = DEXP ( PITSI ) ; VAR2 : = 1.0 / VARl ; ALF : = CST25 • * (VARl - l. IF TSI LSS 1. 30P-15 THEN BEGIN DO ALF := CST25 * ( PITSI » ( 1.0 + 0.5 * PITSI + CST8 * PITSI ) ) ) ( 1.0 END VAR3 VAR4 VAR5 VAR6 VAR7 VAR11 VAR12 VAR13 VAR14 VAR17 VAR18 DEXP(-ALF « ALF ) ; ERF( ALF ) ; 1.0 - VAR4 ; ERF( CST6 * ALF ) ; 1.0 - VAR6 ; ■■ VAR2 * VAR2 ; ■ CST29 « VAR2 ; « 4.0 • CST28 « VAR11 ; : VARl - VAR2 : : -VAR2 + VAR3 ; ■ VAR3 * ( 11.0 * VAR5 - 9.0 ) + 4.0 » VAR5 00105900 00106000 00106100 00106200 00106300 00106400 00106500 00106600 00106700 00106800 00106900 00107000 00107100 00107200 00107300 00107400 00107500 00107600 00107700 00107800 00107900 00108000 00108100 00108200 00108300 00108400 00108500 00108600 00108700 00108800 00108900 00109000 00109100 00109200 00109300 00109400 00109500 00109600 00109700 00109800 00109900 00110000 00110100 00110200 00110300 00110400 00110500 00110600 00110700 00110800 00110900 00111000 00111100 00111200 00111300 00111400 00111500 00111600 00111700 00111800 12 k - 6.0 ; 00111900 VAR19 := VAR3 - 1.0 ; 00112000 VAR21 := 1.0 + VAR11 - 2.0 * VAPb ; 00112100 VAR22 := VAR11 - \/AR3 ; 00112200 VAR23 := VAR5 - VAR3 ; 00112300 IF TSI LSS 1.0PP-15 THEN 00112400 BEGIN 00112500 VAR14 := PITSI » ( 2.0 + CST8 * PITSI • PITSI ) ; 00112600 VAR15 := ALF « ALF ; 00112700 VAR16 := VAR15 « ( - 1.0 ♦ 0.5 » VAR15 » 00112800 (1.0 - CST8 * VAR15 )) ; 00112900 VAR20 := 2.0 » PITSI ; 00113000 VAR24 := VAR20 • ( -1.0 ♦ 0.5 » VAR20 « 00113100 ( 1.0 - CST8 » CST20 )) ; 00113200 VAR21 := VAR24 + 2.0 * VAR4 ; 00113300 IF ALF LSS 1.0SP-7 THEN 00113400 BEGIN 00113500 VAR17 := PITSI * ( 1.0 - 0.5 * PITSI * 00113600 (1.0 + CST8 • PITSI )) + VAR16 ; 00113700 VAR22 := VAR21 - VAR16 ; 00113800 END ; 00113900 END ♦ 00114000 IF ALF LSS 1.000-7 THEN 00114100 BEGIN 00114200 VAR15 := ALF * ALF ; 00114300 VAR16 := VAR15 • ( - 1.0 + 0.5 • VAR15 » 00114400 (1.0 - CST8 * VAR15 )) ; 00114500 VAR18 := 2.0 » VAR16 - VAR4 * ( 15.0 - 11.0 » 00114600 VAR16 ) ; 00114700 VAR19 := VAR16 ; 00114800 VAR23 := - VAR4 - VAR16 ; 00114900 END ; 00115000 LEVEL := VAR14 + CST23 » CST9 * 00115100 ( VAR17 - CST1 * ALF » VAR5 00115200 + CST23 * ( -CST10 * VAR4 + 00115300 CST11 * ALF * ALF * VAR5 - 1.5 * ALF * VAR3 )) 500115400 INJT := ALF * ( CST12 * VAR3 * VAR3 + VAR5 » ( 0.5 » 00115500 VAR4 - CST31) + ALF » ( VAR3 * (-CST14 ♦ VAR5 + 00115600 CST17) + ALF * VAP5 » (CST8 • VAR5 - CST16) ) ) 00115700 + CST2 » ( CST26 » VAR6 + CST15 » 00115800 VAR19 + 0.5 • CST8 * VAP18 ) ; 00115900 VAR8 := 8.0 » CST23 » T * INT I 00116000 VAR9 := VAR5 * ( 3.0 + VAR2 » ( 6.0 • CST23 * ALF 00116100 -2.0 + VAR2 * CST23 » ( ALF « ( 4.0 - 00116200 6.0 * CST23 • ALF ) - CST23 ))) ♦ VAR3 * 00116300 2.0 * CST2 « ( -ALF + CST24 - 2.0 « 00116400 CST23 » VAR2 * ( 1.0 ♦ VAR2 « ( 1.0 - 00116500 1.5 * CST23 * ALF ))) ♦ CST22 * VAR2 * 00116600 VAR2 ? 00116700 VAR10 := 4.0 ♦ T « CST25 * ( VAP5 « ( VAR12 * 00116800 CST13 - VAR13 » CST14 + ALF » ( -6.0 • 00116900 CST16 - VAR13 * CST13 - VAR3 * VAR12 * 00117000 3.0 * CST2 ♦ ALF * ( -3.0 » VAR12 * 00117100 CST16 + VAR3 * ( 2.0 ♦ CST2 ♦ VAR13 « 00117200 CST14 ) + ALF * ( VAP13 « CST16 - 00117300 V/AR5 * VAR13 » CST8 ))) + VAR5 * 00117400 ( -0.5 * VAR12 +ALF * ( 2.0 + 0.5 » 00117500 VAR13 + ALF * VAR12 )) - VAR3 * ( CST2 + 00117600 CST18 » VAR13 )) + VAR3 • ( 2.0 * CST2 * 00117700 125 ( 1.5 + CST30 ) - CST30 • VAR12 - 00117800 VAR13 » CST19 + VAR3 * VAR12 • 2.0 « 00117900 CST3 + ALF * ( 8.0 * CST12 ♦ CST9 » 00118000 VAR12 * ( 1.0 + CST12 ) - VAR3 • ( 2.0 • 00118100 CST3 ♦ CST12 * VAR13 ) - ALF * CST2 * 00118200 ( 2.0 ♦ VAR13 » CST16 ))) + 00118300 VAR13 » CST20 * VAR7 - VAR13 * CST21 ) ; 00118400 ^OR I := 2 STEP 1 UNTIL NXENDM DO 00118500 BEGIN 00118600 PSI(IfJ) := ( LEVEL + VAR8 » CU) ) * Stll ; 00118700 ZETAdvJ] := ( VAR9 + VAR10 » CEI) ) * Stll; 00118800 FDZTII) := PI * ( 8.0 « CST2 * CST29 ♦ 7.0 - 00118900 CST9 • CST25 + REYN « ( CST2 • CST29 * 00119000 ( 1.0 ♦ CST30 ) - 2.0 ) * CIIl ) » 00119100 SU1 ; 00119200 SL0PEII1 := -FDZTU) • CST3 ; 00119300 END ; 00119400 TSI := TSI + DTSI ; 00119500 END ; 00119600 TIME1 := T ; 00119700 END ; 00119800 %»«»««««*«'»«»*«»**»»e«««»#*«#«»#*»»»*«**»»*»*«*»»#*»H*««-»-»'t»»*#*»*##«<»* 00119900 % 00120000 PROCEDURE 3C ; 00120100 BEGIN 00120200 INTEGER L»I»J; 00120300 FOR J := 1 STEP 1 UNTIL NYEND DO 00120400 BEGIN 00120500 ^OR I := 1 STEP 1 UNTIL NXEND DO 00120600 BEGIN 00120700 KOREUtJ] := PSHIfJ] *» 00120800 END ; 00120900 END ; 00121000 END ; 00121100 %tt»«»#»»«»**#»»»*»»««»*»*H»»ttH*tt»»»*«»»»»*»»»»»tt»«»^«>»«tt«H»tt«*»#«#««»«»'tt»00121200 % 00121300 PROCEDURE IMITIALIZENET ; 00121400 % INITIALIZE THE PSI NET WITH THE POTENTIAL FLOW00121500 % SOLUTION* RELAX IT* COMPUTE THE SURFACE 00121600 % VORTICITY*AND SET THE REST OF THE VORTIC I TYOO 121700 * NET EQUAL TO ZERO 00121800 BEGIN 00121900 INTEGER NXP* NXM, II ♦ I»J*K.L ; 00122000 REAL TCG» ETCG, RETCGS» FTG» DIE* SI 1 * PlYlNt CSIGM. 00122100 PIDEL* PIE ; 00122200 LABEL LI ; 00122300 IF SECTION NEO 1 THEN BEGIN FETCHPSIZETA ; BC ; GO TO Li;END; 00122400 NXM := NXEND / 2 ; 00122500 NXP := NXM + 1 ; 00122600 SIGMA := 0.0 ; 00122700 FOR J := 1 STEP 1 UNTIL NYENDM DO 00122800 BEGIN 00122900 TCG := PI • SIGMA ; 00123000 ETCG := EXP(TCG) ; 00123100 RETCGS := 1.0 / (ETCG » ETCG) ? 00123200 FTG := ETCG • (1.0 - RETCGS) ; 00123300 ETA := 0.0 ; 00123400 FOR I := 1 STEP 1 UNTIL NXP DO 00123500 BEGIN 00123600 PIE := PI » ETA ; 00123700 126 SI1 := SIN(PIE) ; PSHItJ] := FTG • SI1 ; ETA := ETA + DELX ; END ; sigma := sigma + dely ; end; piyin := pi * yinf ; csigm := exp(piyin) ; ftg := csigm ; if potentials then BEGIN RETCGS := 1 .0/ (CSIGM*CSIGM) ; F"TG := CSIGM * (1.0 - RETCGS) ; END ; PIDEL := PI * DELX ; FOR I := 1 STEP 1 UNTIL NXP DO BEGIN FI := i-i ; »si(i» nyend] := ftg » sin(pidel » end; for j := 1 step 1 until nyend do BEGIN FOR I := 1 STEP 1 UNTIL NXM DO BEGIN II := NXENDP - I ; PSHII* J] := PSIIIt J] ; END ; end ; for j := 1 step 1 until nyend do BEGIN FOR I := 1 STEP 1 UNTIL NXEND DO BEGIN ZETAt I, J) := 0.0 ; END ; END ; IF PP=1 THEN PRINT1 ( PSI • IPRNT ) ; IF STARTTIME EOL 0.0 THEN BEGIN BC;STREAM; Li: IF STARTTIME GTR 1.0P9-30 AND SECTION EQL 1 THEN SHORTTIMESOLUTION ( END ; * PROCEDURE MATINV < A ,N,S»M, DETERM ) ; DOUBLE ARRAY A[»»»] , B[*»«] ; INTEGER N»M t DOUBLE DETERM ; BEGIN DEFINE IROW = JROW ** ICOLUM = JCOLUM SWAP = T # ; DOUBLE ARRAY PIVOT(0:51 ; DOUBLE T ; INTEGER ARRAY IPI VOT [ : 5 J ♦ INDEX [ : 5t : 2 ] INTEGER JROW» JCOLUM, LI »I»J»K»L ; LABEL HOP720 ♦ HOP777 ? DETERM := 1.0 ; FOR J := 1 STEP 1 UNTIL N DO 3EGIN IPIVOTtJl := ; END ; FOR I := 1 STEP 1 UNTIL N DO FI) surfacevorticity;end; starttime) ; #, AMAX = T «• 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 123800 123900 124000 124100 124200 124300 124400 124500 124600 124700 124800 124900 125000 125100 125200 125300 125400 125500 125600 125700 125800 125900 126000 126100 126200 126300 126400 126500 126600 126700 126800 126900 127000 127100 127200 127300 127400 127500 127600 127700 127800 127900 128000 128100 128200 128300 128400 128500 128600 128700 128800 128900 129000 129100 129200 129300 129400 129500 129600 129700 127 HOP100: HOP105: HOP260: BEGIN LABEL HOP260, HOP380, HOP550 I AMAX := 0.0 ; FOR J := 1 STEP 1 UNTIL N DO BEGIN LABEL HOP105 ; IF (IPIVOTI JJ-1) EQL THEN GO FOR K := 1 STEP 1 UNTIL N 00 BEGIN LABEL HOP100 ; IF (IPIVOTlK] -1) GTR THEN IF (IPIVOTlK] -1) EQL THEN IF (DABS(AMAX) - DABS(AU,K] GO TO HOP100 J ; = k ; A[J,ki ; TO HOP105 ; GO TO H0P777 ; GO TO HOP100 ; ) ) GTR 0.0 THEN ; HOP380: IROW : ICOLUM AMAX : END ; END ; IPIVOTI ICOLUM] := IPIVOTI ICOLUM] ♦ IF (IROW - ICOLUM) EQL THEN GO DETERM := - DETERM I 1 ; TO HOP260 ; FOR STEP 1 UNTIL N DO ; IF FOR L := 1 BEGIN SWAP := AlIROW,L] ; AIIROW,L) := A[ ICOLUM, LJ Al ICOLUM, L) := SWAP ; END ; M LEQ THEN GO TO HOP260 L := 1 STEP 1 UNTIL M DO BEGIN SWAP := BI IROW,L 1 ; BIIROW,L) := BUC0LUM,L1 ; Bt ICOLUM, LI := SWAP ; end; indexii,1j := irow ; indexu,2] := icolum ; pivotii] := a[ icolum, icolum) ; determ := determ « pivotii] IF PIVOTII] EQL AIICOLUM, ICOLUM] c"OR IF FOR FOR STEP THEN := 1.0 ; 1 UNTIL N GO TO DO HOP720 := A[ICOLUM,L] / PIVOTII] GO TO UNTIL HOP380 M DO / pivoti i ) ; L := 1 BEGIN ACIC0LUM,L) END ; . M LEQ THEN L := 1 STEP 1 BEGIN BIIC0LUM,L] := B[ICOLUM,L) END ; LI := 1 STEP 1 UNTIL N BEGIN IF (Ll-ICOLUM) EQL THEN T := AIL1, ICOLUM] ; AIL1, ICOLUM] := 0.0 ; FOR L := 1 STEP 1 UNTIL N DO BEGIN A[L1,L1 := AIL1,L1 - ACICOLUM,L) » T ; END ; IF M LEQ THEN GO TO HOP550 ; DO GO TO HOP550 ; 00129800 00129900 00130000 00130100 00130200 00130300 00130400 00130500 00130600 00130700 00130800 00130900 00131000 00131100 00131200 00131300 00131400 00131500 00131600 00131700 00131800 00131900 00132000 00132100 00132200 00132300 00132^00 00132500 00132600 00132700 00132800 00132900 00133000 00133100 00133200 00133300 00133400 00133500 00133600 00133700 00133800 00133900 00134000 00134100 00134200 00134300 00134400 00134500 00134600 00134700 00134800 00134900 00135000 00135100 00135200 00135300 00135400 00135500 00135600 00135700 128 FOR L := 1 STEP 1 UNTIL M DO BEGIN BIL1.L1 := BIL1»L) - B[ICOLUM,Ll END ; end; HOP550: END ; FOR I := 1 STEP 1 UNTIL N DO BEGIN LABEL HOP710 ; L := N+ 1 - I ; IF (INDEXtL.ll - INDEXtL»21) EQL THEN GO JROW := INDEXtLtl) ; JCOLUM := INDEX[L»2) ; F"OR K := 1 STEP 1 UNTIL N DO BEGIN SWAP := AIK»JR0W] ; A[K*JROW] := A[K*JCOLUM) ; A(K»JCOLUMJ := SWAP ; END ; hop710: go to h0p777; end; HOP720: 3EGIN FORMAT FMT15 ("MATRIX IS SINGULAR" ) ; WRITE (LINE* FMT15 ) ; GO TO hopend ; end ; H0P777: BEGIN END ; END ; % PROCEDURE LSQ1 ( X, Y . W »N*L* A»B»M) ; DOUBLE ARRAY A[***l » B[**»] ; INTEGER N»M*L ; ARRAY X[«lt Y[*»»3* W[»] ; BEGIN DOUBLE ARRAY C(0:5* 0:5] ; DOUBLE DETERM ; INTEGER I* J* K ; FOR I := 1 STEP 1 UNTIL N DO BEGIN ct l*n := l.o ; end ; for j := 2 step 1 until m do BEGIN FOR I := 1 STEP 1 UNTIL N DO BEGIN CtI»J] := CCI.J-II * xt n ; end ; 00135800 00135900 *T ; 00136000 00136100 00136200 00136300 00136400 00136500 00136600 00136700 TO HOP710; 00136800 00136900 00137000 00137100 00137200 00137300 00137400 00137500 00137600 00137700 00137800 00137900 00138000 00138100 00138200 00138300 00138400 00138500 00138600 «««#«»««««««"»«»Q0138700 00138R00 00138900 00139000 00139100 00139200 00139300 00139400 00139500 00139600 00139700 00139800 00139900 00140000 00140100 00140200 00140300 00140400 00140500 00140600 129 end; 00140700 for i := 1 step 1 until m do 00140800 3EGIN 00140900 FOR J := 1 STEP 1 UNTIL M DO 00141000 BEGIN 00141100 A[I»J] := 0.0 ; 00141200 FOR K := 1 STEP 1 UNTIL N DO 00141300 BEGIN 00141400 A(I*J] := All, J] + CCK.I1*C(K,J] ♦ W[K) ; 00141500 END ; 00141600 END ; 00141700 END ; 00141800 FOR I := 1 STEP 1 UNTIL M DO 00141900 3EGIN 00142000 9tltl] := 0.0 ; 00142100 FOR K := 1 STEP 1 UNTIL N DO 00142200 BEGIN 00142300 Bl I» 1 1 := BI I* 1 1 + CtK,I) * YtK»l]*W[Kl ; 00142400 END i 00142500 END ; 00142600 MATING ( AvMfBfLt DETERM ) ; 00142700 END ; 00142800 ^ft«##«»#«ft«««»««tt«tt««««««o«««tt««»4#««»#»«»«»tf#«««tttttt«««««««#»««#«»«««««00 142900 % 00143000 PROCEDURE ORESS ; 00143100 BEGIN 00143200 REAL ARRAY RADI0:5)» WI0:51, AZETA1 : 5* : 11 ♦ ZETASCOtXDl, 00143300 PREStOtXDl ; 00143400 REAL oiSIGt REDEL* THETA, PAREA. CD * SAPEA ; 00143500 INTEGER I* J» K «L ; 00143600 LABEL LI ; 00143700 FORMAT FMT12 ( X50 • "PRESSURE DISTRIBUTION"// X45»"TIME = "♦ 00143800 F"7.3» X6» "NITER ="♦ 15 // X50» "ANGLE'S X9» 00143900 "PRESSURE"/ X50f ••*»»♦»"♦ X9» "««»»»««»" // 00144000 ( X50, F5.1, X10» F7.3) ) ; * 104 00144100 FORMAT FMT13 ( //// X56»"CDP ="» F9.3// X56,"CDF ="♦ 00144200 F9.3 // X56»"CD ="»F9.3) ; * 109 00144300 FORMAT FMT72 ( X50» F5.1, X10* F7.3 ) I 00144400 IF STARTTIME GTR 0.0 AND STARTTIME EOL TIME1 THEN GO TO LU 00144500 SIGMA := 0.0 ; 00144600 FOR J := 1 STEP 1 UNTIL 5 DO 00144700 BEGIN 00144800 °ISIG := PI * SIGMA ; 00144900 WIJ1 := 1.0 ; 00145000 RADCJ1 := EXP(PISIG) ; 00145100 SIGMA := SIGMA + OELY ; 00145200 END ; 00145300 FOR I := 1 STEP 1 UNTIL NXEND DO 00145400 BEGIN 00145500 FOR J := 1 STEP 1 UNTIL 5 DO 00145600 BEGIN 00145700 AZETA[J»1):= ZETAU»J] ; 00145800 END ; 00145900 LSQ1 ( RAD*AZETA«tf*4tltA*Bt3 ) ; 00146000 SLOPEtll := B(2»ll + 2.0 • B[3»1J ; 00146100 END; 00146200 Li: PRESm := 0.0 ; 00146300 REDEL := (4.0/REYN) » (PI* DELX/3.0) ; 00146400 FOR I := 3 STEP 2 UNTIL NXEND DO 00146500 BEGIN 00146600 130 D RESm := PREStI-2] - REDEL * (SL0PEII-21 + 4.0 « 00146700 SLOPEU-11 + SLOPECI) ) ; 00146800 END ; 00146900 FOR I := 1 STEP 2 UNTIL NXEND DO 00147000 BEGIN 00147100 FI := 1-1 ; 00147200 THETA := FI * PI » DELX ; 00147300 ZETAStll := -PRESIII » COS(THETA) ; 00147400 END ; 00147500 PAREA := 0.0 ; 00147600 L := 5 ; 00147700 IF NYENDM MOD 4 NEQ THEN 00147800 BEGIN 00147900 PAREA := 1.50 « (ZETASdl + ZETASC3J) ; 00148000 L := 7 ; 00148100 EMD ; 00148200 FOR I := L STEP 4 UNTIL NXEND DO 00148300 BEGIN 00148400 OAREA := PAREA + (ZETASU-4] + 4.0 » ZETASU-21 00148500 + ZETASU] ) * 00148600 END ; 00148700 PAREA := (PI * 2.0 * DELX/3.0 ) « PAREA ; 00148800 ZETASl U := 0.0 ; 00148900 FOR I := 2 STEP 1 UNTIL NXEND DO 00149000 BEGIN 00149100 fl := 1-1 ; 00149200 THETA:= FI » PI * DELX ; 00149300 ZETASUl := ZETAUvl) • SIN(THETA) ; 00149400 END ; 00149500 SAREA := 0.0 ; 00149600 FOR I := 3 STEP 2 UNTIL NXEND DO 00149700 BEGIN 00149800 SAREA := SAREA + < ZETAStI-2] + 4.0 » ZETASl 1-1 ] 00149900 + ZETASII] ) ; 00150000 END ; 00150100 SAREA := -REDEL « SAREA ; 00150200 CD := PAREA + SAREA ; 00150300 WRITE ( LINE [SKIP 11); 00150400 PGNUM := PGNUM + 1 ; 00150500 WRITE (LINE* FMT75* PGNUM ) ; 00150600 WRITE (LINE* 16 ♦ T ITLEARRAY [ * ] ) ; 00150700 WRITE (LINE* FMT90 ) ; 00150800 WRITE (LINE* FMT12* TIME1»NITER* FOR I := 1 STEP 1 00150900 UNTIL NPRES DO { ANGLE! I ]» PRES(2*I-11 ) ) 5 00151000 IF EVEN = 1 THEN WRITE(LINE* FMT72* ANGLEINPRES + 1] » 00151100 PRES[2»NPRES1 ) ; 00151200 WRITE (LINE* FMT13* PAREA* SAREA* CD ); 00151300 NPP := ; 00151400 END J 00151500 % 00151700 PROCEDURE =RINTPARAMETERS ; 00151800 BEGIN 00151900 PGNUM := PGNUM + 1 ; 00152000 WRITE (LINE* FMT75* PGNUM ) ; 00152100 WRITE (LINE* 16 * TITLEARRAY [ * ] ) ; 00152200 WRITE (LINE* FMT16 ) ; 00152300 WRITE (LINE* FMT4) 5 00152400 % REYN::= REYNOLDS NUMBER 00152500 % YINF::= LARGEST VALUE OF TS1*THE RADIAL CO- 00152600 131 % ORDINATE % LNP::= THE NUMBER OF TIME STEPS BETWEEN PR % OUTS MINUS ONE. % NITER: := TIME STEP COUNTER* INITIALLY ZERO % LNITER::= TOTAL NUMBER OF TIME STEPS FOR T % RUN % TIME2::= MAXIMUM PROCESSOR TIME IN SECONDS % DTIME ::= DIMENSIONLESS TIME INCREMENT WRITE (LINE. % NXEND ::= NUMBER OF MES % CIRCUMFERENTIAL D % NYEND ::= NUMBER OF MES % RADIAL DIRECTION % WFP ::= RELAXATION PARA % TIME1 : := INITIAL DIMEN % POTENTIAL : := IF It OUT * POTENTIAL FLOWt EL FMT5. REYNt YINFt LNPt LNPP DTIMEt TIMElt NITE POTENTIAL. NXEND. IOY.DELXt DELY.STA DSV.PPZ.IRECt SECTION ) ; H POINTS IN THE IRECTION H POINTS IN THE METER SIONLESS TIME ER BC TAKEN AS SE PARALLEL FLOW t LNITERt TIME2. R. MYEND.IQX. RTTIME.SAVE. END ; % PROCEDURE BEGIN IF PZ IF PP IF OS NP : = END ; % PROCEDURE BEGIN IF NP IF SA LOCK( LOCK( END ; DRINTPSIZETA ; OUTPUTRESULTS ; =1 THEN PRINT1 ( ZETA. ZPRNT ) ; =1 THEN PRINT1 ( PSI. PPRNT ) ; V EQL 1 THEN SAVEPSIZETA ; ; finalprintroutine ; neq then begin outputresults; press end ; ve=i then begin;storepsizeta;storeconstants;end; SAVEP) ; SAVEZ) i % PROCEDURE BEGIN LABEL SETUP GETPS OUTPU PRESS H0P165tNITER NP : = NPP : TIME1 GETPS H0P165. HOP190 ; GET ; IZETA ; TRESULTS ; := NITER + 1 ; np + l ; = npp + l ; := TIMEl + DTIME ; IZETA ; * PRINT CURRENT VORTICITY AND STREAMFUNCT ION VALUES * AND COMPUTE THE PRESSURE DISTRIBUTION % WHENEVER REQUESTED BY LNP IF NP EQL LNP THEN OUTPUTRESULTS; 00 INT 00 00 00 HIS 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 *»««*«»tt»»»#*«»»tt«<»«»00 % PROCEDURE MAIN ; BEGIN LABEL H0P165* HOP190 ; INITIALIZEPROGRAM ; COMPUTECONSTANTS ; PRINTPARAMETERS ; IF PPZ EQL 1 THEN BEGIN PRINTPSIZETA ; GO TO HOP190;END; INITI4LIZENET ; IF SECTION NFQ 1 THEN GO TO H0P165 ; OUTPUTRESULTS : PRESS ; % MAIN PROGRAM LOOP H0P165:NITER := NITER + 1 ; NP := NP + 1 ; NPP := NPP + 1 ; TIME1 := TIME1 + DTIME ; % COMPUTE THE NEW VORTICITY BY THE % ALTERNATING DIRECTION IMPLICIT METHOD VORTE ; VORTS ; % RELAX THE STREAM FUNCTION WITH THE % ACCELERATED LIEBMAN METHOD STREAM ; % COMPUTE THE SURFACE VORTICITY SURFACEVORTICITY ; * PRINT CURRENT VORTICITY AND STREAMFUNCTION VALUES * AND COMPUTE THE PRESSURE DISTRIBUTION % WHENEVER REQUESTED BY LNP IF NP EQL LNP THEN OUTPUTRESULTS ; IF NP 3 EQL LNPP THEN PRESS ; % IF PROGRAM RUNS OVERTIME GO TO END OF JOB CHECKTIMELIMIT ; % IF DESIRED NUMBER OF TIME STEPS HAVE % BEEN COMPLETED THEN GO TO END OF JOB IF NITER GEQ LNITER THEN GO TO HOP190 ; 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 % OTHERWISE COMPUTE THE NEXT TIME STEP GO TO H0P165 ; hopi9o:end ; % 00 main ; 00 H0P999:FINAL°RINTR0UTINE ; 00 HOPENDtTYMOUT ; 00 END ; 00 %H>^«»*«O-»«««4»»»«*«tt»«H*«»»»tt*««»^*«'tt«tt«««tt»*»»»4»«»»tt»«»»«tttttt«tt»»«-tt«-«#i*««'00 % 00 158700 158800 158900 159000 159100 159200 159300 159400 159500 159600 159700 159800 159900 160000 160100 160200 160300 160400 160500 160600 160700 160800 160900 161000 161100 161200 161300 161400 161500 161600 161700 161800 161900 162000 162100 162200 162300 162400 162500 162600 162700 162800 162900 163000 163100 163200 163300 163400 163500 163600 163700 163800 163900 164000 164100 164200 164300 164400 164500 164600 END. 133 00164700 13^ APPENDIX E CICHY/UNS BEGIN DECLARATION OF GLOBAL VARIABLES RE REAL SUBROUTINE ENTE^SYSTEMPAWAMETERS : 3EGIN MREY := 40.0 .; TLMT := O.OB : LIMIT := 0.00001 : NR1 MX2 : = MX3 : = \IX4 : = MX5 : = ME1 = NY? = \JY3 = NJY4 = MY5 = =>INC t — i ; EMD ! R£ REAL SL MPOUTINK INITIALIZECOORUINATES ; 3EGU J r lLL Rl WITH 2.0. 1.9, 1.8, 1.7. 1 .b. 1.5. 1.4, 1.3. 1.22, 1.160, 1.11637, 1.0774. ) .0515* 1 .03424. 1 .0226, 1.01497. L. 00985. 1 .00643. 1.00416. 1.00264. 1.00162S. l.o ; WITH 1 .O009S. 1 .0005. 1.0002. r lLL X? -11.0. -10.0. -9.0. -8.0. -7.0* -6.0. -S.0. -4.24b. -3.69?, -3. 2MB. -2.99?, -2.775. -2.616, -2.5. -2.4. -2.3* -2.2* -2.1. -2.0. -1.97'3JB. -1.90211* -1,78201. -1.61803. -1.41421. -1 .17557, -0.9079H. -0.61803. -0.31287. 0.0. +0.31287 5 r lLL X3 wITH -0.31?87, 0.0. ♦0.312B7. 0.<*44 IS AT X = ^.4 X4[0) : = X3I135) : X4(0) IS AT X=9.9 135 % % % % % % LOOP j:=l, 1,151 DO X4U51) IS X4[J] := X4IJ-1] AT X = 25 X5C0] := X4[ 149] ; X5(0) IS AT X=24.8 LOOP j:«l*l»126 DO X5IJ] := X5IJ-11 X5U26) IS AT X=50 + 0.1 + 0.2 REGION 1 HAS REGION 2 HAS REGION 3 HAS REGION 4 HAS 25 MESH ROWS* Rl= 2 TO Rl= 1 29 MESH ROWS* X2=-ll TO X2= 136 MESH KOWS, X3= Tu X3= 10 151 MESH ROWS* X4 = 10 TU X4= 25 REGION 5 HAS 126 MESH ROWS* X5= 25 TO X5= 50 ^ILL EK01 r lLL EH 1 ] FILL Y2I01 r lLL Y310] WITH 0.0* +1.95* 0.2* 0.45* 0.b2406* 0.70483* 0.76414. 0.81463* 0.86027* 0.90301* 0.94401. 0.98408* WITH 0.0* 1.0* 02390. 06409* 10536* 14858* 19500* 1.24682* 1.30882* 1.39892* 1.6* 1.85* WITH +14.0* 10.0* 6.0* 3.288* 2.5* 2.1* 1.78201* 0.90798* -0.31287* -1.41421* -1.97538* -2.3. -2.775* -4.245* -8.0* -12.0* WITH 0.0* 0.0* 0.0* 0.25* 0.5* 0.64357* 0.71775, 0.77477* 0.82407* 0.86901, 0.91132* 0.V5207* 0.99204. 0.0* 1.00796* 1.03188* 1.07224* 1 .11381* 1.15755* 1.20483* 1.25816* 1.32340* 1.42867* 1.65* 1.9* 13.0, 9.0* 5.0* 2.992* 2.4, 2.0, 1.61803, 0.61803, -0.61803, -1.61803, -2.0, -2.4, -2.992. -5.0, -9.0, •13.0, 0.0, 0.0, 0.05, 0.3, 0.55, 0.66086, 0.73005, 0.78510, 0.83333, 0.87764, 0.91957, 0.96011, 1.0 ; 0, 01592, 03989, 08043, 12236, 16667, 1.21490, 1 .26995, 33914, 450, 7, 95, 12.0, 8.0, 4.245, 2.775, 2.3, 1.97538, 1.41421, 0.31287, -0.90798, -1.78201, -2.1, -2.5, -3.288, -6.U, ■10.0, ■14.0 ; 0.0, 0, 1, 35, 57133 67660 74184 79517 84245 ,88618 0.92776 0.96812 0.0, 0.15, 0.4, 0.60108, 0.69118, 0.75318, 0.80500, 0.85142, 0.89464, 0.93591, 0.97610, 0.0, 1.04793, 1 .08868, 1.130*9, 1.17593, 1.22523, 1.28225, 1.356^3, 1.50, 1.75, 0.0, 11.0* 7.0* 3.69^* 2.616* 1.90211* 1.17557* 0.0* -1.17557* -1.90211* -2.616* -3.69^, -7.0* ■11.0* 0.0* 1.05600* 1.09699, 1.13973, 1.18537, 1.23586. 1.29517, 1 .37594, 1.55, 1 .8, 0.05 : 0.0 0.0 136 ^ILL F ILL r lLL r lLL +14.0, 13.0, 12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.6225, 4.245, 3.9685, 3.692, 3.490, 3.288, 3.140, 2.992, 2.8835, 2.775, 2.6955, 2.616, 2.556 2.5, 2.45, 2.4, 2.35, 2.3, 2.25, 2.2, 2.15, 2.1, 2.05, 2.0, 1.97538, 1.95, 1.90, 1.85, 1.60, 1.75, 1.70, 1 .65, 1.60, 1.55, 1.50, 1.45, 1.40. 1.35, 1.30, 1.25 ; Y3[ 1 ] WITH 0.0, 0.0, 0.0, 0.0, 0.0, 1.20, 1.15, 1.10, 1.05, 1.00, 0.95, 0.90, 0.65, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05, 0.0, -0.05, -0.10, -0.15, -0.20, -0.25, -0.30, -0.35, -0.40, -0.45, -0.50, -0.55, -0.60, -0.65, -0.70, -0.75, -0.80, -0.65, -0.90, -0.95, -1.00, -1.05, -1.10. -1. 15, -1.20, -1.25 ; Y3t2J WITH 0.0, 0.0, 0.0, 0.0, 0.0, -1.30, -1.35, -1.40, -1.45, -1.50. -1.55, -1.60, -1.65, -1.70 ; -1.75. -1.80, -1.85, -1.90, -1.95, -1.97538, -2.00, -2.05, -2.10, -2.15. -2.20. -2.25, -2.30, -2.35. -2.40, -2.45, -2.50, -2.556, -2.61*, -2.6955, -2.775, -2.8835, -2.9Sy2, -3.140, -3.288, -3.490, -3.692, -3.966b, -4.245, -4.6225, -5.0, -6.0, -7.0, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0 ; Y4[01 WITH 0.0, 0.0, 0.0, 0.0, 0.0, +14.0, 13.0, 12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.6225, 4.245, 3. 9b65, 3.692, 3.490, 3.288, 3.140, 3.00, 2.90, 2.8, 2. /, 2.6, 2.5, 2.4, 2.3, 2.2, 2.1, 2.0, 1.9, 1.8, 1.7, 1.6, 1.5, 1.4, 1.3, 1.2, 1.1, 1.0, 0.9, 0.8, 0.7, 0.6, 0.5, o.<+, 0.3, 0.2, 0.1, o.o ; Y4[l] WITH 0.0, O.U, 0.0, 0.0, 0.0, -0.1, -0.2, -0.3, -0.4, -0.5, -0.6 -0. U -O.o, -0.9, -1 .0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.6, -1.9, -2.0, -2.1, -?,iio -2.3, -2.4, -2.5, -2.6, -2.7, -2.6, -2.9, -3.0 -3.140, -3.288, -3.490, -3.692, -3.9685, -4.245, -4.8225, -5.0, -6.0, -7.0, -8.0, 137 -9.0» -14.0 ; TILL Y5fO] WITH +14.0, 10.0, 5.0» 3.288. 2.2, 1.2, 0.2, -0.8, -1.8, -3.9685, -8.0, -13.0, -10.0, 13.0, 9.0, 4.245, 3.0, 2.0, 1.0, 0.0, -1.0, -2.0, -3.0, -4.245, -9.0, -14.0 -11.0, 12.0, 8.0, 3.9685, 2.8, 1.8, 0.8, ■0.2, •1 .2, •3.288, ■5.0, ■10.0, -12.0, 11.0, 7.0, 3.692, 2.6, 1.6, 0.6, ■0.4, ■1.4, ■2.4, ■3.490, ■6.0, ■11. 0, -13. U, 6.0, 3.4*0, 2.4, 1.4, O.'f, ■0.6, ■1 .6, •2.b, •3.692, ■7.0, ■12.0, COMPJTECONST 3EGIN DEFINE A PE REAL SUBROUTINE AnTS ; R 6. + DEFINE AL & + R & IN P I HE S.9 EN DEFINE Rt beg: LOOf DEFINE H REG LOO r lLL r ILL r lLL FILL r lLL r lLL FILL r ILL ^ILL ~ILL "ILL 'ILL r lLL r lLL MHO] AM1[ 1 ] AM2[01 81 82 8 3 84 85 Alt AH A2( A3! A3t A3t A4[ A<+[ A5[ L IN P IL BE E'Ni rn MIJ AM Th IT IS VA FI Wl TH WITH WITH WI Th WITH 0] W w w ■H w M «l W = AP = AW = Ak Al ^ 6.AH Al & 5»Alf II 6, l:=& GIN MIIL ** II 6, :=M GIN P[ IL I.) u a E CE ST h ONG AT T IS SOL LUES LL S I a: = 1 .0/ (S.A1I 6.IA] ) ) nu ; ia: = i.o/<6,ah MAI) ) un ; LIM \BM 6,81: = 5.IA1 » (RTR<1,,6.A1[6.IA)) UA] * (RTL (l,,6,Ai[6.IA)) II,1,6.L1M 00 ) := 1 . 0/ C&dl [ IL 3 # (J.HK IL-1 1+6,611 IL1) ) ; ? LIM KUP 6.81: = I»l, ; AND CB ; (bbbSSSbbS N(OAAAAAAA ANtrtOOOOOO E AN (800000 EAM ( OOUliOO ROL CONSTANTS FOR EACH ROW LUOE ALL °OINTS INCIUL CIRCLt FFFFFFFFFFFFF0(16) ) ; oooood6)) ; E PATTERNS FOP THE PART ERE THE CYLINDRICAL OT NCB ; ) OP CB ; SRb5S5S(16) ) : AAAAAAAAA ( 16) ) ; ooooooooo?< 16) ) OOUOOOuOUO (16) ) ooooooonu?) PE KEAL SUBROUTINE PE REAL SUBROUTI MEWPSI (PCPOINT PSI, PCPOI PCPOINT BM» PCPO PINT IY» PINT IX CI NT S ) ; BEGIN MODE :- PMOOE ; TEST := PSICJ] : ID := 1.0 / (BMC IX ] + NE NT ZETA, PIimT J, PINT K, INT HP, PCPOINT AM, PCPOINT AP » CI NT OUT LIM, BOOLEAN PMODE* *PC IX] + AMC IY1 + APC IY] ) ; 139 PSTfJ] := ID * (O.B » 7ETA[J] + RMtlX] » RTR(lf»PSI[Kl) + HPIIX] * RTLU»»PSIIK]) + AM[IY1 * PSIIJ-SJ + AP(IY) <* PSHJ + S] ) : OIF := PSKJ) - TEbT ; IF ABS(OIF) GTR LIMIT THEN LIM := 1 ; END ; % THE ITERATIVE SOLUTION OF THE STREAM FUNCTION % EQUATION! CYCLtS THROUGH THE 4 SECTIONS PE REAL SUPROUTINE STREAM1 ; BEGIN % TREAT THE INFLOW BOUNDARY MODE := M?fO) ; PSAfO] := psa[ i i ; % SET UP THE TOP A NO BOTTOM BOUNDARIES MODE := M2BCL ; LOOP J.* = l»l»28 DO BEGIN PSAtJJ := RTL( 1»»PSA( J]) + 1.0 ; END ; MODE := M?BCR : LOOP j:=l»l»28 DO BEGIN PSAU) := RTR(1»»PSAC JJ) + 1.0 : END : % % SWEEP THE ODD POINTS IN THE MESH Swl := CEVEN? ; SW? := CODD? pass := l ; A := l ; PASS?: S := 1 : * NUMBER OF pe ROWS D ER MESH ROw IY := o ; LOOP C:=1»?*2H DO BEGIN J := C + SW1 ; * ROWS BEING UPDATEu K := C + SW? ; % ADJACENT POINTS * HOTH RELATIVE TO BUFFER ORIUlN pmouE := M2IC+A1 : IX := j ; NE**PSI (PSA. ZSA* J, K. BM?, HP?. AM?, AP?, IY. IX, LIM, PMOUE. S ) ; END ; IF ^ASS=1 THEN BEGIN PASS := ? : A := ; SW1 := CODD? : Sw? := CEVEN? ; GO TO ^ASS? ; END : % % NOW MOVE PSA[?8] TO A n/ORK AREA, EXPAND IT % BY INTERPOLATION SO IT CORRESPONDS TO THE MESH % ROW WIDTH IN REGION 3. AND INSERT IT IN PSa[3h] END ; PE REAL SUBROUTINE STREAM? ; REG IN END : 140 PE REAL SUBROUTINE STREAM3 ; BEGIN END ; PE REAL SUBROUTINE STREAMS ; BEGIN END : MITER := ; ITERATE: BEGIN LIM := : STREAM1 STPEAM2 STREAM3 STREAM4 IF LIM = THEN GO TO CONVERGED ? NITER := NITER + 1 ; GO TO ITtRATE ; converged: end ; pe real subroutine surfacevorticity ; BEGIN end ; PE PEAL SUBROUTINE INITIALIZESTPEAM ; BEGIN OEFINE INS f>X2 S.Y2 &I2 kJ2; = &Y2UJ2] * (1.0 - 1.0 / START X = 9.9 ENU X=11.B BEGIN LOOP w:=0.1.1 DO BEGIN MODE := M4tW] ; PSBtK + W] := INS X4 Y4 N w : ZSBtK+wi := 0.0 ; END ; N := n + l ; END ; % WRITE SECTION 2 ON DISK AND % READ SECTION <♦ IN TO BUFFER t) END ; PE REAL SUBROUTINE INITIAL3 ; BEGIN N := N-2 ; LOOP K:=0. 2.243 DO * START X=11.7 END X=23.H BEGIN LOOP w: =0.1.1 UO BEGIN MODE := M41W] » PSA [K + w] := INS X4 Y4 N W ; ZSA [K+W] := 0.0 : END ; N := N + 1 ; END ; % WRITE SECTION 3 ON DISK AND % READ SECTION 1 INTO BUFFER A END ; PE REAL SUBROUTINE INITIALS ; Ite BEGIN N := N - 2 ; LOOP K:=0»?»27 DO % START X=23.7 END X=25 BEGIN LOOP W:=0»1»1 DO BEGIN MODE := M4[W] ; PSBfK+W] := INS X4 Y4 N W ; ZS9IK+W] := o.o ; END ; N := N + 1 ; END ; M := o ; LOOP K:=30»l»155 DO % START X=24.8 ENDS X=49.b BEGIN MOOE := ms ; PSBCK] := INS XS Yb N ; ZSBtK] := o.o ; N := n + l ; END ; N := ; LOOP K: = 1 START » = a ENDS R=i BEGIN LOOP w:=0»l»l DO BEGIN MODE := MllW] ; PSBfK+W] := RUN) « SIN(F1[W1) : ,