Main Subroutines of FASTEST/CALCULATION

Subroutine celnbuvw

Description/Structure:

This routine is represents the discretization of the convective and diffusive fluxes of the momentum equations for the cell faces on block boundaries. For that reason it is quite similar to the routine celuvw. Again, the implementation is designed such that the routine is able to compute the corresponding coefficients for the east, north and top face of a control volume. However, the following documentation refers to the east face only.

The routine starts with the storage of some variables of minor importance. The following loop is designed such that all computational node having an west face on the block boundary are addressed. Then, the identification numbers of the neighboring nodes are determined by employing the corresponding offsets (ine, inn...). Additionally, the interpolation factors for CDS and TBI are stored (fxe, fxw, tfp...).

Now, the computation of the diffusive flux starts. For this purpose, three vectors are determined. The first vector (dxks, dyks, dzks) points from the the geometrical location of node P to its neighboring node E. The second (dxet, dyet, dzet) is a connection of the midpoints of the south and north edge of the east face. Finally, the third (dxzd, dyzd, dzzd) represents the corresponding connection for the top and bottom edge. The vector product of the second and third vector yields a vector representing the cell face area (b11, b12, b13).

Then, several coefficients are computed which arise from the transformation of the derivatives of the velocities in the diffusive terms. In principle, they represents the Jacobian of the applied transformation (vole, bt11, bt12, bt13...) and the finite difference approximation of the derivatives (duns, dutb, dvns...). Then the part of the diffusive term which is treated implicitly is stored in game. The remainder of the diffusive fluxes is then summed up in sueu, suev, suew for the $u$-, $v$- and $w$-equation, respectively. The splitting is done for two reasons. First, the diffusive flux contains terms which can not be treated implicitly. Second, the splitting allows to use the same coefficient matrix for all momentum equations. Finally, the implicit part of the diffusive flux is multiplied with the face area.

Now, the convective flux is computed. Note, that only the first order part of the discretization is treated implicitly. The CDS part goes into the source term. To perform the first order discretization (Upwind), the direction of mass flux is determined such that ce or cw contain the absolute value of the mass flux. ae1 or aw1 represent parts of the explicitly treated part of the convective fluxes.

Finally, the matrix entry is computed as a sum of the convective and diffusive terms. Since the east face of a control volume is also the west face of the neighboring volume, the computed fluxes hold also for that computational node. Thus, the corresponding matrix entry can be determined by considering the corresponding terms. Note that due to the first order discretization of the convective terms both matrix entries are not identical.

Eventually, the explicit parts of convection and diffusion are summed up in the corresponding source terms for the $u$-, $v$- and $w$-equations (su, sv, sw) if the block is connected to another block, i.e. the corresponding cell face is an internal face. Otherwise, the explicit part of the diffusion is added to the source terms in the case of wall, outlet or symmetry boundary condition. To be honest, I have no idea why.

Subroutine celp2

Description/Structure:

In principle, this subroutine represents the discretization of the pressure correction equation. The implementation is designed such that the routine is able to compute the corresponding coefficients for the east, north and top face of a control volume. However, the following documentation refers to the east face only.

The routine starts with setting a logical variable such it contains the information if the computations are carried for a coarse grid level or not. Then , a loop over all control volumes having an internal east face is initialized. Then, the identification numbers of the neighbouring nodes are determined by employing the corresponding offsets (ine, inn, ...). Additionally, the interplation factors for CDS and TBI are stored (fxe, fxw, tfp, ...).

Next, two vectors are computed which connect the the midpoints of the south/north edge and the top/bottom edge, respectively. The vector product of both yields the a vector representing the cell face area (arex, ary, arz). The the interpolated value of the inverse of the central coefficient of the momentum equations (apu, apv, apw) multiplied by the corresponding component of the cell face area vector is stored in due, dve, dwe. Additionally, the density is interpolated and a pressure difference is computed. Then, the velocities ue, ve, we are computed at the midpoint of the east face by aemplying the pressure weighted interpolation. On a coarse grid level (multi grid), only the correction of the velocity components are considered such that the starting value of the velocities u1, v1 ,w1 is interpolated and substracted. Instead of employing the coarse grid velocities directly, the restricted fine grid mass flux is employed (fcfad).

Next, the coefficents of the system matrix are computed. Since the east face of a control volume is also the west face of the neighbouring volume, the computed values holds also for the west coefficient of the neighbouring computational node. Additionally, the mass flux over the east face is computed.

Subroutine caluvw

Description/Structure:

The subroutine caluvw, together with calcp forms the basis of the SIMPLE-Algorithm in FASTEST. In caluvw the momentum equations are assembled and solved for $u$, $v$, and $w$. The following discussion will (at least for now) concentrate on the core functionality and omit ``add-ons'' like e.g. porous media or turbulence modelling, since they are not essential for the understanding of caluvw's workings. In the course of the subroutine, the coefficients of the equations are first initialized and then the contributions of various effects are added onto them. Finally, the resulting system of equations is solved.

In the first part of caluvw a number of vectors and variables are initialized. These are (in order of appearance) the vectors $b[ewtbnsp]$, which are auxiliary vectors used for various purposes, and res, in which the residuals from the SIP-solver are stored. The exception in this part of the code is the calculation of the pressure gradients across the CV (calcdp). Initialization continues with the underrelaxation factors (variables starting with urf). Next, all per-block initializations are done. This happens at the beginning of a loop over all blocks and starts with the source terms $su$, $sv$, and $sw$. Depending on the value of loldfl (which is true only in the Crank-Nicolson-Scheme in prepstep.F) the source terms are initialized. Lastly, the central coefficients $ap$, $spv$, and $sp$ are set to zero.

Now that the coefficients have sensible initial values, various effects (buoyancy, etc.) are taken into account by incrementing the coefficients accordingly. The main contribution to the matrix coefficients are due to the convective and diffusive fluxes. These fluxes are both calculated and added to the coefficients in the routines celuvw (c.f. p. ) and celnbuvw (c.f. p. ).

Next, by calling moduvw, boundary conditions are applied.

As a last step before the solver is called, the components of the unsteady terms are added to the matrix coefficients and the source term vectors.

Now that all coefficients have their final values, the resulting systems of equations are solved for $u$, $v$, and $w$.

Subroutine celuvw

Description/Structure:

The subroutine celuvw represents the discretization of the convective and diffusive fluxes of the momentum equations. The routine is designed so that it can compute the corresponding coefficients for the east, north and top face of a control volume. In the following however, we refer to the east face only. The north and top face follow by analogy. Due to a special treatment of the diffusive fluxes, the computed matrix holds for the $u$-, $v$- and $w$-equation.

The routine starts with storing the flux blending coefficient as gam. Next, a loop over all control volumes with an internal east face begins. Then, the indices of the neighboring nodes are assigned to ine, inn, etc. and the interpolation factors for CDS and TBI are stored in fxe, fxw, tf[pent].

Now, the computation of the diffusive fluxes starts. For this purpose, three vectors are determined. The first vector (dxks, dyks, dzks) points from the the geometrical location of node P to its neighbouring node E, the second (dxet, dyet, dzet) is a connection of the midpoints of the south and north edge of the east face, and the third (dxzd, dyzd, dzzd) represents the corresponding connection for the top and bottom edge. The vector product of the second and third vector yields a vector representing the cell face area (b11, b12, b13).

Several auxiliary coefficients are computed, which arise from the transformation of the velocity derivatives in the diffusive terms. Basically they represent the Jacobian of the applied transformation (vole, bt11, bt12, bt13,...) and the finite difference approximation of the derivatives (duns, dutb, dvns,...). The nature of the discretized diffusive fluxes poses two problems: First, they contain terms which cannot be treated implicitly without creating additional non-zero diagonals in the system matrices. Second, their fully implicit handling leads to different matrices for $u$, $v$, and $w$. Thus, the diffusive fluxes are split into an explicit and an implicit part. The latter is stored in game, the former is summed up in sueu, suev, suew for the $u$-, $v$- and $w$-equation, respectively. Concluding the calculation of the implicit part of the diffusive flux, it is multiplied by the cell face area.

Now the convective flux is computed. Note that only its first order UDS discretization is treated implicitly. The direction of the mass flux is taken into account by assigning its absolute value to either ce or cw. To achieve second order accuracy, the difference between the UDS and the CDS discretizations is stored in ae1 and aw1 and later added to the source terms (explicit treatment).

(

Finally, the matrix entry is computed as a sum of the convective and diffusive terms.

Since the east face of a control volume is also the west face of the neighbouring volume, the computed fluxes hold also for that comutational node.

Thus, the corresponding matrix entry can be determined by considering the corresponding terms.

Note that due to the first order discretization of the convective terms both matrix entries are not identical.

)

Finally, the explicit parts of convection and diffusion are summed up in the corresponding source terms for the $u$-, $v$- and $w$-equation (su, sv, sw).

Subroutine fmg3d

fmg3d (full multigrid 3d) is the real main routine of FASTEST, it controls the execution of the FASTEST code.

The basic procedure is as following:

• Definition of the ring and farm communication structures required for data input (initco1(), p. ).

• Input of file names and file handling (fnames(), p. )

• Initialize error handling (initerr(), p. )

• Check status of some files and delete them if necessary and supported (filestat(), p. ).

• Input of control values (readinp(), p. )

• Read input file for visualization (readvis(), p. )

• Read file for start values for all fluid regions, default file: id/PRJNAME.prop (readprop(), p. )

• Input of grid datas (readmap(), p. )

• Data input and initializations (modinp(), p. )

• Initialize values for the SIP-solver (setsip(), p. )

• Define index vector necessary for vectorization (setvec(), p. 

• Open some output files (openfls(), p. )

• Preparing time steps.

• Start of full multigrid loop.

• Interpolation of the solution to a next finer grid, if full multigrid is employed (prepgrid(), p. )

• Time loop:
  • Output of cfl-, diffusion- and cell-Reynolds number (cfldif(), p. )

  • V-cycle loop:
    • Start loop over finest grid on actual muligrid-level

    • Calculate inflow flows and residuum normalization factors (flxin(), p. )

    • Calculate momentum equations (caluvw(), p. ).

    • Calculate outlet boundary conditions for velocities and correct fluxes (outbcvel(), p. ).

    • Actualization of variable fluid properties (iprops(), p. )

    • Calculate pressure equations, incl. boundary conditions, etc. (calcp(), p. )

    • Correct $u$, $v$ and $w$ on corner lines (corvel(), p. )

    • Check of convergence and output of residuals and monitoring values.

    • END of loop over finest multigridlevel -- if convergence criteria are satisfied: jump to the next multigrid level

    • Set boundary pressure (extrapolates pressure from interior points) (bpres(), p. ) (to the block-faces or to the boundaries???)

    • Restriction and solution from finest to coarsest grid (restr()), and prolongation and solution from coarsest to finest grid (inter()).

  • Actualization of pressure boundary values (bpres(), p. ).

  • Output of data for postprocessing or continuation of a computation and of information of the solution process (computing time, numbers of iterations,...). (wwrite(), p. ).

• Postprocessing -> (wwrite(), p. , outsurf(), p. , surfdat(), p. , fornumeca(), p. , fortec(), p. ,...

• END of Full-Multigrid-loop

• Close files and save data for restart

A short description of the solution algorithm is in sec. .