HPC MSU
    Cover
    Preface
  1. Introduction
  2. Boundary-Conforming Coordinate Systems
  3. Transformation Relations
  4. Numerical Implementation
  5. Truncation Error
  6. Elliptic Generation Systems
  7. Parabolic and Hyperbolic Generation Systems
  8. Algebraic Generation Systems
  9. Orthogonal Systems
  10. Conformal Mapping
  11. Adaptive Grids
  12. Appendix A
    Appendix B
    Appendix C
    1. Code Development Exercises
    2. Computer Exercises
    3. Listing of Routine for Computer Exercises
    4. Examples for Computer Exercises
    References

    Downloadable Version (PDF)
Numerical Grid Generation
Foundations and Applications
By: Joe E. Thompson, Z.U.A. Warsi and C. Wayne Mastin

APPENDIX C

CODE DEVELOPMENT AND COMPUTER EXERCISES

1. Code Development Exercises

1. Take the computational region to be a rectangle on which the curvilinear coordinates are defined to be equal to the indices of the field arrays, i.e., = I and = J, with x = X (I,J) and y = Y (I,J).

Make provision for reading in values of x and y on any segments of the boundary of the computational region. Generate x and y in the interior by interpolating linearly between the top and bottom boundaries. Plot the grid.

2. Modify the code to allow horizontal (in the computational plane) interpolation, as well, the choice being specified by input.

3. Now add the choice of interpolation from the four corners (tensor product interpolation).

4. Finally add the choice of transfinite interpolation.

5. Generalize the interpolation to cubic Hermite interpolation, with the grid being orthogonal at the boundary.

6. Generalize the interpolation to use the hyperbolic tangent distribution function, rather than being linear, with specified relative spacing on each end.

7. Modify the code to provide for reading in x and y on any segment of any horizontal or vertical line in the computational region. Also provide for the interpolation tion to be done on any rectangular segment of the computational region (including a segment that is only a line.)

8. Add another field array TYPE (I,J) which is a flag to identify each point as one for which the x,y values are (1) fixed, e.g., specified points on the physical boundary, (2) out of the computation, e.g., points inside a slab, or (3) to be generated. Provide for the designation (1) and (2) to be made by input for any rectangular segment of the computational region, the default being to the designation (3).

9. Modify the dimensions of the field arrays so that an extra layer of points surrounds the computational region. Also add two more field arrays, ILINK (I,J) and JLINK (I,J). Provide for any segments of any horizontal or vertical lines to be designated as image points in TYPE by input, i.e., points for which the values of x and y are set equal to those at some other point. Also provide for the indices of these other points to be put in ILINK and JLINK by input.

10. Add an elliptic generator, based on Laplace equation, to the code. Use the algebraic generator (the interpolation) to provide the initial guess for point SOR iteration.

11. Add control functions to the elliptic generator. Let the control function be evaluated on the boundaries and interpolated into the field by transfinite interpolation.

2. Computer Exercises

1. Generate an algebraic grid between two concentric circles. Use linear interpolation between the circles.

2. Generate an algebraic grid between two ellipses, both of which are centered at the origin but which may have different eccentricities, using interpolation between the ellipses. Compare grids generated using linear and Hermite interpolation, the latter being orthogonal at the boundaries.

3. Generate a C-type algebraic grid for an ellipse inside an outer boundary formed by a semicircle replacing one side of a rectangle:

Compare (1) vertical interpolation in the computational region boundary, (2) horizontal interpolation, (3) tensor product interpolation, and (4) transfinite interpolation, using linear interpolation in each case. Note that (2) and (3) are totally unreasonable.

4. Generate an algebraic grid for a circular simply-connected region by (1) unidirectional interpolation, (2) tensor product interpolation, and (3) transfinite interpolation. Note that here only (3)gives a reasonable grid. Compare linear and Hermite interpolation for (3).

5. Repeat Exercise 4 with a triangular boundary.

6. Using the boundary configuration of Exercise 3, but with a hyperbolic tangent point distribution on the right-hand boundary of the physical region with smaller spacing at the centerline than at the top and bottom. Compare algebraic grids generated using (1) linear interpolation between the inner and outer boundaries, (2) nonlinear interpolation, based on the hyperbolic tangent, between the inner and outer boundaries, (3) transfinite interpolation with linear blending functions, and (4) transfinite interpolation using the boundary point distribution (in terms of relative arc length) as the blending functions. Note that only (2) and (4) preserve the boundary point distribution in the field.

7. Generate an algebraic grid for a square inside a rectangle using linear interpolation between the inner and outer boundaries. Note the propagation of the boundary slope discontinuities into the field. Generate a grid from an elliptic generation system for the same boundary point distribution and note the difference.

8. Generate an algebraic grid for a square inside a circle using linear interpolation between the inner and outer boundaries. Show that it is possible to position the points on the circle such that the grid overlaps the corners of the square. Generate a grid from an elliptic generation system for the same boundary point distribution and note the difference.

3. Listing of Routine for Computer Exercises

      SUBROUTINE INTERP 
      PARAMETER(NI=20,NJ=20,N=5) 
      COMMON/COORD/X(NI,NJ),Y(NI,NJ) 
      COMMON/CONST/CHOICE,IMAX,JMAX,NA,DS1,DS2 
      COMMON/ATTR/IAL(N),IAX(N),IAY(N),JAL(M),JAX(N),JAY(N) 
      COMMON/COEF/AI(N),BI(N),CI(N),DI(N),AJ(N),BJ(N) 
      COMMON/COEF/CJ(N),DJ(N) 
      DIMENSION P(NI,NJ),Q(NI,NJ),XX(0:NI,O:NJ),YY(0:MI,O:NJ) 
      DIMENSION X1(NI,NJ),X2(NI,NJ),Y1(NI,NJ),Y2(NI,NJ) 
      INTEGER CHOICE 
C  
C     BOUNDARY INTERPOLATION  
C  
C      X          X ARRAY OF XI-ETA COORDINATE  
C      Y          Y ARRAY OF XI-ETA COORDINATE  
C      IMAX       MAX. NUMBER OF GRID IN XI AXIS  
C      JMAX       MAX. NUMBER OF GRID IN ETA AXIS  
C      NA         MAX. NUMBER OF ATTRACTIONS  
C      DSl        SPECIFIED LENGTH OF INITIAL INTERVAL  
C      DS2        SPECIFIED LENGTH OF FINAL INTERVAL   
C      ATTR       ARRAY OF ATTRACTION TO LINES/POINTS  
C      COEF       ARRAY OF COEFFICIENT FOR ATTRACTION  
C  
C      CHOICE    
C         1       VERTICAL INTERPOLATION     
C         2       HORIZONTAL INTERPOLATION   
C         3       TENSOR PRODUCT INTERPOLATION   
C         4       TRANSFINITE INTERPOLATION      
C         5       HERMITE CUBIC INTERPOLATION    
C         6       HYPERBOLIC TANGENT INTERPOLATION    
C         7       ELLIPTIC GRID GENERATION ( SOR ITERATION )       
C         8       ATTRACTION TO COORDINATES           
C             
       IF(CHOICE.EQ.1) GO TO 100  
       IF(CHOICE.EQ.2) GO TO 200   
       IF(CHOICE.EQ.3) GO TO 300                      
       IF(CHOICE.EQ.H) GO TO 400 
       IF(CHOICE.EQ.5) GO TO 500                            
       IF(CHOICE.EQ.6) GO TO 600  
       IF(CHOICE.EQ.7) GO TO 700                                  
       IF(CHOICE.EQ.S) GO TO 800  
C  
C      **** VERTICAL INTERPOLATION ****   
C             
 100   DO 110 I 1,IMAX        
       DO 110 J 1,JMAX     
       RJ1=FLOAT(JMAX-J)/FLOAT(JMAX-1)  
       RJ2=FLOAT(J-1)/FLOAT(JMAX-1)     
C      *** ( EQ. 8-1 )  
       X(I,J)=RJ1*X(I,l)+RJ2*X(I,JMAX)            
 110   Y(I,J) RJl*Y(I.1)+RJ2*Y(I,JMAX)         
       RETURN                                     
C                                                  
C      **** HORIZONTAL INTERPOLATION ****             
C                                                  
 200   DO 210 I=l,JMAX                              
       DO 210 J=1,IMAX                            
       RI1=FLOAT(IMAX-I)/FLOAT(IMAX-1)            
       RI2=FLOAT(I-1)/FLOAT(IMAX-1)   
C      *** ( EQ. 8-1 )                                
       X(I,J)=RI1*X(1,J)+RI2*X(IMAX,J)            
 210   Y(I,J)=RI1*Y(1,J)+RI2*Y(IMAX,J)              
       RETURN                                     
C                                                  
C      **** TENSOR PRODUCT INTERPOLATION ****         
C                                                   
 300   DO 310 I=1,IMAX                              
       DO 310 J=l,JMAX                            
       RI1=FLOAT(IMAX-I)/FLOAT(IMAX-1)            
       RI2=FLOAT(I-1)/FLOAT(IMAX-1)               
       RJl=FLOAT(JMAX-J)/FLOAT(JMAX-1)            
       RJ2=FLOAT(J-1)/FLOAT(JMAX-1)               
C      *** ( EQ. 8-69 )                               
       X(I,J)=RI1*RJ1*X(1,1)+RI1*RJ2*X(1,JMAX) 
      *+RI2*RJ1*X(IMAX,1)+RI2*RJ2*X(IMAX,JMAX) 
       Y(I,J)=RI1*RJ1*Y(1,1)+RI1*RJ2*Y(1,JMAX) 
      *+RI2*RJ1*Y(IMAX,1)+RI2*RJ2*Y(IMAX,JMAX) 
 310   CONTINUE                                      
       RETURN                                      
C 
C      **** TRANSFINITE INTERPOLATION ****             
C 
 400   DO 410 I=l,IMAX                               
       DO 410 J=1,JMAX                             
       RI1=FLOAT(I-1)/FLOAT(IMAX-1)                
       RI2=FLOAT(IMAX-I)/FLOAT(IMAX-1)             
       X1(I,J)=RI1*X(IMAX,J)+RI2*X(1,J)            
 410   Yl(I,J)=RI1*Y(IMAX,J)+RI2*Y(1,J)              
       DO 420 I=1,IMAX 
       DO 420 J=1,JMAX 
       RJ1=FLOAT(J-1)/FLOAT(JMAX-1) 
       RJ2=FLOAT(JMAX-J)/FLOAT(JMAX-1) 
       X2(I,J)=RJ1*(X(I,JMAX)-X1(I,JMAX))+RJ2*(X(I,1)-X1(I,1)) 
 420   Y2(I,J)=RJ1*(Y(I,JMAX)-Y1(I,JMAX))+RJ2*(Y(I,l)-Y1(I,1)) 
C      *** ( EQ. 8-73 ) 
       DO 430 I=1,IMAX                                         
       DO 430 J=1,JMAX                                         
       X(I,J)=X1(I,J)+X2(I,J)                                  
 430   Y(I,J)=Y1(I,J)+Y2(I,J)    
       IF(CHOICE.NE.4) GO TO 740                               
       RETURN                                                  
C                                                               
C      ***  HERMITE CUBIC INTERPOLATION (ORTHOGONAL BOUNDARY) *** 
C                                                               
 500   DO 510 I=1,IMAX                                          
       DO 510 J=1,JMAX                                     
       XX(I,J)=X(I,J)                                 
 510   YY(I,J)=Y(I,J)                       
       DO 520 J=1,JMAX                      
       XX(O,J)=XX(IMAX-1,J)            
       YY(O,J)=YY(INAX-1,J)       
       XX(IMAX+1,J)=XX(2,J)  
 520   YY(IMAX+1,J)=YY(2,J)  
       DO 530 I=1,IMAX       
       DO 530 J=1,JMAX       
       RJJ=FLOAT(J-1)/FLOAT(JMAX-1)
C      *** ( EQ. 8-6 a and b, n=2 )   
       PHI1=(1.+2.*RJJ)*(1.-RJJ)*(1.-RJJ) 
       PHI2=(3.-2.*RJJ)*RJJ*RJJ           
       PSI1=(1.-RJJ)*(1.-RJJ)*RJJ         
       PSI2=(RJJ-1.)*RJJ*RJJ              
C  
C      ** CAL. NORMAL DERIV. **                
C  
       XXI1=.5*(XX(I+1,1)-XX(I-1,1))      
       XXI2=.5*(XX(I+1,JMAX)-XX(I-1,JMAX))
       YXI1=.5*(YY(I+1,1)-YY(I-1,1))      
       YXI2=.5*(YY(I+1,JMAX)-YY(I-1,JMAX))
       UNIT1=SQRT(XXI1*XXI1+YXI1*YXI1)    
       UNIT2=SQRT(XXI2*XXI2+YXI2*YXI2)  
C      *** ( EQ. 3-108 )                    
       XNl=-YXI1/UNIT1*DSl           
       XN2=-YXI2/UNIT2*DS2      
       YN1=XXI1/UNIT1*DSl  
       YN2=XXI2/UNIT2*DS2  
C      *** ( EQ. 8-5 )       
       XX(I,J)=PHI1*XX(I,l)+PHI2*XX(I,JMAX)+PSI1*XN1+PSI2*XN2 
 530   YY(I,J)=PHI1*YY(I,1)+PHI2*YY(I,JMAX)+PSI1*YN1+PSI2*YN2 
       DO 540 I=1,IMAX  
       DO 540 J=1,JMAX 
       X(I,J)=XX(I,J) 
 540   Y(I,J)=YY(I,J)   
       RETURN                                                
C 
C      **** HYPERBOLIC TANGENT SPACING INTERPOLATION **** 
C   
 600   TOL=1.0E-10    
C      *** ( EQ. 8-49, 50 and 51 
       A=SQRT(DS2/DS1)                                       
       B=1./(FLOAT(JMAX-1)*SQRT(DS1*DS2))                    
C      *** INITIAL GUESS BY SERIES EXPANSION                    
       DELTA=SQRT(6.*(B-1.))                
       DO 610 IT=1,20                    
       RESID=SINH(DELTA)/(DELTA*B)-1.         
       IF(ABS(RESID)LT.TOL) GO TO 630      
 610   CALL AITKEN(DELTA,RESID,DELTO,RO,RSO)  
       PRINT 620, RESID,DELTA,IT-1  
 620   FORMAT(//, 5X, 'DELTA IS NOT CONVERGE ?', 5X, 2E15.5, 
      *5X, I3, //)                   
       GO TO 660                   
 630   CONTINUE                     
C      *** ( EQ. 8-52, 53 and 54 )    
       DO 650 I=1,IMAX             
       DO 650 J=2,JMAX-1           
       RATIO FLOAT(J-1)/FLOAT(JMAX-1)  
       U=.5*(1.+TANH(DELTA"(RATIO-.5))/TANH(.5*DELTA)) 
       S=U/(A+(1.-A)*U)     
       X(I,J)=X(I,1)+(X(I,JMAX)-X(I,1))*S   
 650   Y(I,J)=X(I,1)+(Y(I,JMAX)-Y(I,1))*S    
 660   RETURN    
C   
C      **** ELLIPTIC GRID GENERATION ( SOR ITERATION ) ****  
C      *** CAL. P AND Q ON THE BOUNDARY **   
C   
 700   DO 710 I=1,IMAX    
       DO 710 J=1,JMAX   
       XXI=.5*(X(I+1,J)-X(I-1,J))       
       XXIXI=X(I+1,J)-2.*X(I,J)+X(I-1,J)  
       XETA=.5*(X(I,J+1)-X(I,J-1))    
       XETA2=X(I,J+1)-2.*X(I,J)+X(I,J-1)  
       YXI=.5*(X(I+1,J)-Y(I-1,J))         
       XXIXI=Y(I+1,J)-2.*X(I,J)+Y(I-1,J)  
       YETA=.5*(Y(I,J+1)-Y(I,J-1))        
       YETA2-Y(I,J+1)-2.*Y(I,J)+Y(I,J-1)  
       IF(ABS(XETA2).LT.10E-3) XETA2=0.   
       IF(ABS(YETA2).LT.10E-3) YETA2=0.   
       RXI2=XXI*XXI+YXI*YXI  
       RETA2=XETA*XETA+YETA*YETA                               
C      *** ( EQ. 8-70 )   
       P(I,J)=(XXI*XXIXI+XXI*YXIXI)/RXI2     
       Q(I,J)=(XETA*XETA2+YETA*YETA2)/RETA2  
 710   CONTINUE   
C  
C      **    INTERPOLATE P AND Q BETWEEN BOUNDARY **  
C            P : VERTICAL, Q : HORIZONTAL        
C  
       DO 720 I=1,IMAX       
       DO 720 J=l,JMAX          
       RJ1=FLOAT(JMAX-J)/FLOAT(JMAX-1)             
       RJ2=FLOAT(J-1)/FLOAT(JMAX-1)                   
       RI1=FLOAT(IMAX-I)/FLOAT(IMAX-1)                   
       RI2=FLOAT(I-1)/FLOAT(IMAX-1)                         
       P(I,J)=RJl*P(I,l)+RJ2*P(I,JMAX)                         
       Q(I,J)=RI1*Q(1,J)+RI2*Q(IMAX,J)                         
 720   CONTINUE 
C  
C      **    INITIAL GUESS WITH TRANSFINITE INTERPOLATION **  
C      
       GO TO 400  
 740   CONTINUE  
C  
C      *** ITERATION ( SOR ) ***  
C  
       ITMAX=200 
       TOL=10.E-5                                              
       W=1.8                                                   
       DO 760 IT=1,ITMAX                                       
       ERRX=0.                                                 
       ERRY=0.                                                 
       DO 750 J=2,JMAX-1                                       
       DO 750 I=2,IMAX-1                                       
       XXI=.5*(X(I+1,J)-X(I-1,J))                              
       YXI=.5*(Y(I+1,J)-Y(I-1,J))                      
       XXIXI=X(I+1,J)+X(I-1,J)                                 
       YXIXI=Y(I+1,J)+Y(I-1,J)                                 
       XETA=.5*(X(I,J+1)-X(I,J-1))                             
       YETA .5*(Y(I,J+1)-Y(I,J-1))                             
       XXIETA=.25*(X(I+1,J+1)-X(I+1,J-1)-X(I-1,J+1)+X(I-1,J-l))
       YXIETA=.25*(Y(I+1,J+1)-Y(I+1,J-1)-Y(I-1,J+1)+Y(I-1,J-1))
       XETA2=X(I,J+1)+X(I,J-1)                                  
       YETA2=Y(I,J+1)+Y(I,J-1)                                  
C      *** ( EQ. 6-18 and 6-20 )    
       G11=XXI*XXI+YXI*YXI                                      
       G22=XETA*XETA+YETA*YETA 
       G12=XXI*XETA+YXI*YETA 
       XTEMP=.5*(G22*(P(I,J)*XXI+XXIXI)+G11*(Q(I,J)*XETA+XETA2) 
      *-2.*G12*XXIETA)/(G11+G22)                                     
       YTENP=.5*(G22*(P(I,J)*YXI+YXIXI)+G11*(Q(I,J)*YETA+YETA2) 
      *-2.*Gf2*YXIETA)/(Gll+G22)                                    
       XTENP=W*XTEMP+(1.-W)*X(I,J)                                  
       YTEMP=W*YTEMP+(1.-W)*Y(I,J)                              
       ERRX=AMAXO(ERRX,ABS(XTEMP-X(I,J)))                           
       ERRY=AMAXO(ERRY,ABS(YTEMP-Y(I,J)))                           
       X(I,J)=XTEMP                                                 
       Y(I,J)=YTEMP                                                 
 750   CONTINUE            
       IF(ERRX.LT.TOL.AND.ERRY.LT.TOL) GO TO 78O                    
 760   CONTINUE                 
       PRINT 770,ERRX,ERRY,IT-1                                     
 770   FORMAT(//, 5X, 'X AND Y ARE NOT CONVERGE ?', 2E15.5,         
      *5X, I5, //)                                                  
 780   CONTINUE 
       IF(CHOICE.EQ.8) GO TO 830                                    
       RETURN   
C 
C      **** ATTRACTION TO COORDINATE LINE/POINT **** 
C 
 800   DO 810 I=1,IMAX  
       DO 810 J=1,JMAX                                              
       P(I,J)=0.                                            
 810   Q(I,J)=0.    
       DO 820 NS=1,NA                                      
       DO 820 I=1,IMAX                                  
       DO 820 J=1,JMAX                               
       XL=FLOAT(I-IAL(NS))                            
       XI=FLOAT(I-IAX(NS))                         
       XJ=FLOAT(J-IAY(NS))                      
       YL=FLOAT(J-JAL(NS))                   
       YI=FLOAT(I-JAX(NS))               
       YJ=FLOAT(J-JAY(NS))            
C      *** ( EQ. 6-30 )         
       P(I,J)=P(I,J)-AI(NS)*(XL/ABS(XL))*EXP(-CI(NS)*ABS(XL))  
      *-BI(NS)*(XI/ABS(XI))*EXP(-DI(NS)*SQRT(XI*XI+XJ*XJ))     
       Q(I,J)=Q(I,J)-AJ(NS)*(YL/ABS(YL))*EXP(-CJ(NS)*ABS(YL))  
      *-BJ(NS)*(YJ/ABS(YJ))*EXP(-DJ(NS)*SQHT(YI*YI+YJ*YJ))     
 820   CONTINUE  
       GO TO 400 
 830   CONTINUE  
       RETURN  
       END  

4. Examples for Computer Exercises