Code development¶

The current project is for developing elliptic grid generator in 2-dimensional domain. Hereafter, the program developed in this project is called ‘GridGen’.

GridGen Code summary¶

The present project is to make a grid-generator for 2-D computational domain around a modified NACA 00xx series airfoil in a channel.

The source code contains the following directories:

io - input/output related routines
main - main program driver
math - thomas algorithm
modules - main grid generator solver routines
utils - list of useful FORTRAN utilities used within the program
gridpy - python wrapper for gridgen main program

Also a ‘CMakeLists.txt’ file is also included for cmake compiling.

$ cd GridGen/src/
$ ls
$ CMakeLists.txt  io  main math modules utils gridpy

The io folder has io.F90 file which contains ReadInput(inputData) subroutine. It also includes input_file_xml which describes the structure of the user run-time input file located in the main ‘src’ directory, and output.F90 for storing data in bothb Tecplot and Python format.

The main folder is only used for containing the code driver file. The main routines is run by gridgen.F90 which calls important subroutines from the rest of folders.

Details of GridGen development¶

The GridGen code is made for creating 2-D computational domain with pre-described points value along the 2D airfoil geometry. The schematic below shows the flow chart of how the GridGen code runs.

The source code shown below is gridgen.F90 and it calls skeletal subroutines for generating grid structure. The main features of the main code is to (1) read input file, (2) make initialized variable arrays, (3) set initial algebraic grid points, (4) create elliptic grid points, and (5) finally write output files:

PROGRAM main

      USE xml_data_input_file
      USE GridSetup_m,ONLY:InitGrid,GridInternal,EndVars
      USE parameters_m,ONLY:wp
      USE SimVars_m,ONLY:fileLength
      USE GridTransform_m,ONLY: InitArrays,CalcCoeff,PhiPsi,TriDiag,Jacobian,&
              EndArrays
      USE output_m,ONLY:WritePlotFile,WriteRMS

      IMPLICIT NONE

      TYPE(input_type_t) :: inputData
      CHARACTER(LEN=fileLength) ::output= 'none'
      CHARACTER(LEN=fileLength) :: rmsout = 'rmslog.dat'
      INTEGER :: i
      REAL(KIND=wp) :: rms

      CALL InitGrid(inputData)
      CALL InitArrays
      CALL Jacobian
      CALL WritePlotFile('InitGrid','"x","y","Jacobian","Phi","Psi"',inputData)
      IF (inputData%setup%iControl == 1) CALL PhiPsi
      WRITE(*,*) 'Solving Elliptical Grid...............'
      DO i=1,inputData%setup%nmax
                WRITE(*,*) 'IT = ',i
                CALL CalcCoeff
                CALL TriDiag(rms)
                CALL WriteRMS(i,rms,rmsout)
                IF (rms <= inputData%setup%RMSres) EXIT
      ENDDO
      CALL Jacobian
      CALL WritePlotFile(output,'"x","y","Jacobian","Phi","Psi"',inputData)
      CALL EndArrays
      CALL EndVars

END PROGRAM main

Creation of algebraic grid points¶

The code starts to run by reading the important input parameters defined by the user in the input_file.xml file. The input data file first contains all the details geometrical description of our grid points. Then the code reads airfoil geometry data from this input file, which provides the bottom edge points of the domain. The input file also contains four vertex points in \((x,y)\) coordinates. Thus those points forms a 2-dimensional surface. Based on these boundary grid points, the code runs with Algebratic grid generating subroutine and gives initial conditions for elliptic solution for grid transformation.

The gridgen.F90 file first refers to InitGrid subroutine defined in GridSetup.F90 file. The main function of this routine is to call again multiple subroutines defined in same file. The subroutine definition shown below summarizes the how the code runs for the grid initialization:

SUBROUTINE InitGrid(inputData)

        USE xml_data_input_file
        USE io_m, ONLY: ReadInput

        IMPLICIT NONE

        TYPE(input_type_t) :: inputData

        CALL ReadInput(inputData)
        CALL InitVars()
        CALL BottomEdge()
        CALL SetBCs()
        CALL GridInternal()

END SUBROUTINE InitGrid

ReadGridInput: Reads user defined variables and parameters for grid configuration.
InitVars: Initialize the single- and multi-dimensional arrays and set their size based on the input parameters.
BottomEdge: Generate point values for airfoil geometry.
SetBcs: Generate grid points along 4 edges of the computational domain.
GridInternal: Based on grid points along the edges and surfaces, this routine will create interior grid points that are aligned with user-defined grid point interpolations.

Creaction of elliptic grid points¶

In order to determine the elliptic grid points with the pre-specified boundary points, the following Poisson equations, which is given in previous Project description section, have to be resolved numerically. The coefficients of the equations can be determined by:

\[ \begin{align}\begin{aligned}{\alpha}=x_{\eta}^{2} + y_{\eta}^{2}\\{\beta}=x_{\xi}x_{\eta} + y_{\xi}y_{\eta}\\{\gamma}=x_{\xi}^{2} + y_{\xi}^{2}\end{aligned}\end{align} \]

Then, applying finite difference approximation to the governing equations can be transformed into the linear system of equations. The arranged matrix form of equations shown below can be solved for unknown implicitly at every pseudo-time level. At every time loop, the code updates the coefficients composed of \(\phi\) and \(\psi\), and adjacent points. The detailed relations of each coefficients are not shown here for brevity.

\[ \begin{align}\begin{aligned}a_{i,j} x_{i-1,j}^{n+1} + b_{i,j} x_{i,j}^{n+1} + c_{i,j} x_{i+1,j}^{n+1} = d_{i,j}\\e_{i,j} y_{i-1,j}^{n+1} + f_{i,j} y_{i,j}^{n+1} + g_{i,j} y_{i+1,j}^{n+1} = h_{i,j}\end{aligned}\end{align} \]

Above equations can be numerically evaluated by the following descritized expressions:

\[ \begin{align}\begin{aligned}a_{i,j} = e_{i,j} = {\alpha}_{\text{ }i,j}^{n} \left(1 - \frac{\phi_{i,j}^{n}}{2} \right)\\b_{i,j} = f_{i,j} = -2 \left({\alpha}_{\text{ }i,j} + {\gamma}_{\text{ }i,j} \right)\\c_{i,j} = g_{i,j} = {\alpha}_{\text{ }i,j}^{n} \left(1 + \frac{\phi_{i,j}^{n}}{2} \right)\\e_{i,j} = \frac{{\beta}_{\text{ }i,j}^{n}}{2} \left(x_{i+1,j}^{n} - x_{i+1,j-1}^{n+1} - x_{i-1,j+1}^{n} - x_{i-1,j-1}^{n+1} \right) - {\gamma}_{\text{ }i,j}^{n} \left( x_{i,j+1}^{n} + x_{i,j-1}^{n+1} \right) - \frac{{\beta}_{\text{ }i,j}^{n}}{2} \psi_{i,j}^{n} \left( x_{i,j+1}^{n} - x_{i,j-1}^{n+1} \right)\\h_{i,j} = \frac{{\beta}_{\text{ }i,j}^{n}}{2} \left(y_{i+1,j}^{n} - y_{i+1,j-1}^{n+1} - y_{i-1,j+1}^{n} - y_{i-1,j-1}^{n+1} \right) - {\gamma}_{\text{ }i,j}^{n} \left( y_{i,j+1}^{n} + y_{i,j-1}^{n+1} \right) - \frac{{\beta}_{\text{ }i,j}^{n}}{2} \psi_{i,j}^{n} \left( y_{i,j+1}^{n} - y_{i,j-1}^{n+1} \right)\end{aligned}\end{align} \]

where \(n\) and \(n+1\) indicate pseudo time index. Thus above equations will update grid point coordinates for \(n+1\) time level by referring to already resolved \(n\) time level solution. Note that the pseudo time looping goes along the successive \(j\)-constant lines. Therefore, when writing the code, time level index in above equations was not considered as a separate program variable because \(j-1\) constant line is already updated in the previous loop.

The expressions above are only evaluted in the interior grid points. The points on the boundaries are evaluated seprately by applying given solutions as problem handout.

Once initial algebraic grid points are created, the code is ready to make elliptic grid points with some control terms in terms of \(\phi\) and \(\psi\). gridgen.F90 file contains the necessary subroutine calls to evaluate the elipptical grid as shown below:

IF (inputData%setup%iControl == 1) CALL PhiPsi
WRITE(*,*) 'Solving Elliptical Grid...............'
DO i=1,inputData%setup%nmax
          WRITE(*,*) 'IT = ',i
          CALL CalcCoeff
          CALL TriDiag(rms)
          CALL WriteRMS(i,rms,rmsout)
          IF (rms <= inputData%setup%RMSres) EXIT
ENDDO

Before going into the main loop for solving poisson equations, the code calculate control terms with \(\phi\) and \(\psi\) only if the user defines the option in the input_file.xml file. Even though the assigned project made an assumption of linear interpolated distribution of \(\phi\) and \(\psi\) at interior points, the GridGen code is designed to allow \(\phi\) and \(\psi\) be weighted in \(j\) and \(i\) directions, respectively. This effect is made by the grid stretching formula.

Here, main DO-loop routine goes with setup of coefficients of governing equations and Thomas loop. The Thomas loop operates with line Gauss-Siedel method for resolving unknown variables, \(x\) and \(y\), with tri-diagonal matrix of coefficients of finite difference approximation equation in a \(j\) = constant line.

RMS residual¶

In order to avoid infinite time-looping for the Thomas method, the GridGen code employs the following definition of RMS residual based on the new (\(n+1\)) and old(\(n\)) values of grid point coordinates.

\[\text{RMS}^{n} = \sqrt{\frac{1}{N} \sum_{i=2}^{imax-1} \sum_{jmax-1}^{j=2} \left[\left(x_{i,j}^{n+1} - x_{i,j}^{n} \right)^{2} + \left(y_{i,j}^{n+1} - y_{i,j}^{n} \right)^{2} \right]}\]

where \(N = 2x(\text{imax}-2) x (\text{jmax}-2)\) and the RMS criterion is user-specified as as small number or by default it is set as: \(1\text{x}10^{-6}\). In this code, the convergend is assumed to be achieved when RMS residual is less than the RMS criterion.

Table Of Contents

Previous topic

Next topic

This Page