- Introduction
- Access to Q1-set information
- Access to flow variables (SOLVEd and STOREd)
- Example of access to flow variables
- Use of the FN subroutines
- Coding for the example
- Direct access to the F-array
- Coding for the example
- Comparison of FN subroutines and direct f-array access
- Spare working space
- Setting fluid properties
- Setting boundary conditions and sources
- Access to geometrical quantities
- Service functions and subroutines
- Summary

The structure of the user-accessible FORTRAN-coding in EARTH is described in a previous lecture in this Course.

In particular, the structure of GROUND was discussed. GROUND is the empty Ground-Station module in which users insert their own coding.

This lecture explains how the insertion of coding can be effected, and the tools provided by PHOENICS to facilitate the task.

Some examples of usage will be given, but plenty of these are also available in the exemplary-Ground subroutine, GREX.

Full details are given in the Encyclopedia, under the entries: FORTRAN Coding, F-array, Function library and Function subroutines.

Each PIL variable has a counterpart (with the same name) in GROUND. These GROUND variables are automatically assigned the Q1-set value.

For instance, the following variables are available for use in GROUND, and have the values set in the Q1 file:

- Logicals:CARTES, ONEPHS, PARAB, BFC, USEGRD, ...
- Integers:NX, NY, NZ, LSWEEP, IZPRF, NXPRIN, ...
- Reals:XULAST, RHO1, RHO2, RHO1A, ENUL, PRESSO, ...
- Character:NAMGRD, NSAVE, ...
- Real Arrays:PRNDTL, PRT, FIINIT, VARMIN, ...

A full list of these variables can be obtained by inspection of the COMMONs in the directory d_includ.

All these variables can be used anywhere in GROUND, and not just in the Group they 'belong' to.

However, the user will not find GROUND equivalents of Q1-variables locally defined by him/herself.

E.g., the following Q1 setting

REAL(THIS); THIS = 3.0

will not imply the availability, in GROUND, of THIS.

Of course, one can always declare the variable (if this is needed) and assign the value in FORTRAN in GROUND. This is bad practice, as it requires recompilation and relinking every time the value is changed.

However, PHOENICS does provide a means of sending values of non-PIL variables from the Q1 file to GROUND.

The Q1-To-Ground Arrays

__Four arrays are provided, both in the PIL (Q1) and FORTRAN (GROUND) for the
transmission of data between both. These are:__

Array Name |
Type |
Default dimension |
PARAMETER |

RG | real | 200 | NRG |

IG | Integer | 200 | NIG |

LG | Logical | 100 | NLG |

CG | Character*4 | 100 | NCG |

In the example of the previous panel, setting in the Q1 file RG(1)=THIS would allow the use of THIS (with the name RG(1)) in GROUND.

Alternatively, putting THIS = RG(1) in Group 1, Section 1 would allow the use of the same variable name as in Q1.

Of course, GROUND would be quite useless without the possibility of
manipulating flow variables (such as velocities or enthalpies) or variable
fluid-properties (such as density or temperature). How this is done in PHOENICS is the
subject of the next panels.

Notes on terminology:

- The term "PIL variables" has been used above to refer to those variables which PHOENICS provides for use in the Q1 file, and which are also available in the FORTRAN coding of GROUND (such as RHO1, LSWEEP, NX).
- The term "flow variables" (rather than just "variables") will be
used from now on to refer to those variables which are SOLVEd for or STOREd through PIL
commands (such as enthalpy or velocity components).

The Flow-Variable Index

The flow variables are identified, in PHOENICS, through an integer index.

In PIL, the following 50 indices are provided to refer to the flow variables:

- For Phase 1: P1, U1, V1, W1, R1, C1, C3, C5, .... C35
- For phase 2: P2, U2, V2, W2, R2, C2, C4, C6, .... C34
- Other:
- RS (shadow-phase volume-fraction),
- KE (kinetic energy) and
- EP (rate of dissipation of KE)

Note however that the assignment of the Ci variables to one phase or the other can be changed through the TERMS command.

These same indices are also available in GROUND as integer variables.

In GROUND, these PIL-variables (P1, U1, etc.) are NOT 2-D or 3-D arrays which hold the field values of each flow-variable.

They are, rather, integer INDICES which are used by PHOENICS to determine WHERE the variable is stored in the central array F.

They can be thought of as POINTERS into the F-array. Their values are always fixed.

Using The Flow-Variable Indices

__PHOENICS provides two ways of retrieving and setting values of the flow
variables, namely:__

- The FN library
- Direct access to the F-array

The first can access slabs of flow-variable values; the second has more general abilities, but it is normally used to access slabs.

These methods will be discussed next using an example.

We are solving, in a single-phase problem, for the z-velocity, W1 and for the
mass fraction of a chemical species, C1. We want to obtain, at each point, the mass flux
of C1.

This is given, at each cell, by the product RHO1 *W1*C1 (the density RHO1 will be assumed to be 1.0).

The mass flux is to be stored in a spare concentration variable, for which the
following provision has been made in the Q1 file:

STORE(MASS)

STORE will search backward for the first unused variable. It will then arrange for storage and printout, changing the default NAME to the one specified.

Where is the coding to be inserted?

The first decision to be made is in which Group/Section of GROUND the coding is
to be inserted.

This depends of course, on the purpose of the coding. The lecture on the structure of Ground gave some indications in this respect.

In this case, the mass flux needs to be computed only once, when all the other variables have their final values. Also, it is only required by the user, not by EARTH. The coding will therefore go into GROUP 19.

Furthermore, the section in which the coding is going to be inserted must be within the
IZ loop, if we want to have slabwise access to the flow variables.

DO ISTEP=1, LSTEP ! Transience loop | Visit Sec 1: start of time-step | DO ISWEEP=1, LSWEEP ! Outer-iteration loop | | Visit Sec 2: start of sweep | | DO IZ=1, NZ ! Single-sweep loop | | | Visit Sec 3: start of IZ slab | | | DO ITHYD=1, LITHYD ! Hydrodynamic loop | | | | Visit Sec 4: start of iteration | | | | Solve equations | | | | Visit Sec 5: finish of iteration | | | END DO | | | Visit Sec 6: finish of IZ slab | | END DO | | Visit Sec 7: finish of sweep | END DO | Visit of Sec 8: finish of time step END DO

Group 19, Section 6 (Finish of IZ slab) is a good candidate for this example, since it is within the IZ loop, as required. Unfortunately it is also within the ISWEEP loop.

Insertion of the coding in Group 19 Section 6 will entail its execution at the end of each sweep, and not only at the end of the last one. This causes no harm, of course; but it is a waste of computing time (but not a great deal).

This could be avoided by placing the coding inside a BLOCK-IF:

IF(ISWEEP.EQ.LSWEEP) THEN

coding

ENDIF

Use Of The FN Library Of Subroutines

PHOENICS provides a set of subroutines which perform various kinds of arithmetic operation with slab values directly on the F-array.

These subroutines are generally referred to as the FN subroutines, and are fully documented in the on-line help-dictionary (entry FUNCT). Just a few examples are listed below.

Nomenclature:

- Y, X, X1, .... represent integer indices of flow-variables (e.g., W1) Y is always the index of the OUTPUT X, X1 are always the pointers of the INPUT Arithmetic operations on Y and X are NOT allowed, as this would result in an invalid index.

- A, B, C, .... represent scalar constants (e.g., 3.14, RHO1A) Arithmetic IS allowed on the scalars.

Examples Of FN Subroutines

__CALL FN1(Y, A)__

will set the current-slab values of the flow-variable Y to A; in short Y=A.

Example: CALL FN1(C2, 0.0)

CALL FN0(Y,X)

will set the current-slab values of the flow-variable Y to those of the flow-variable X; in short, Y=X.

Examples:

- CALL FN0(C2, C1) ! [C2 = C1]
- CALL FN21(Y, X1, X2, A, B) ! [Y = A + B * X1 * X2
- CALL FN27(Y, X) ! [Y = Y / X]
- CALL FN40(Y) ! [Y = ABS(Y)]

The use of the FN routines depends on knowledge of the index values. In the
example case, we know the indices for velocity (W1) and concentration (C1), but we do not
know which index was chosen by STORE for the variable we called MASS.

The integer function LBNAME returns the index associated with a given name. Thus LBNAME('MASS') will give us the missing index.

Hence, the required coding, in Group 19 Section 6 is:

196 CONTINUE

C * ----------SECTION 6---------- FINISH OF IZ SLAB

C FN21(Y, X1, X2, A,B) Y=A+B*X1*X2

CALL FN21 (LBNAME('MASS'), W1, C1, 0.0, 1.0)

RETURN

If the Q1 had said NAME(C2)=MASS; STORE(MASS), we could have used C2.

Concatenated FN Calls

Of course, concatenated calls to FN subroutines can be performed for complex
operations for which an FN subroutine does not exist. The example above could have been
coded as:

C.....FN0(Y,X).......Y=X

CALL FN0(C2, C1)

C.....FN26(Y, X).....Y=Y*X

CALL FN26(C2, W1)

PHOENICS users can also create their own FN subroutines. This is recommended if the same sequences of FN calls appears many times in a GROUND routine.

Sample listings of FN routines, which can be used as templates, are provided in the Encyclopedia, under FNxxx.

The F-array is the [large] central array in which most of the information is kept. Particularly, the cell values of the flow variables are stored in the F-array.

PHOENICS provides a tool-kit of functions which facilitate access to the values stored in F.

The F-array and its tool-kit are fully described in the Encyclopedia, under F-ARRAY. Only a sub-set of that information will be presented in this lecture.

How The Information Is Stored In F

The cell-values of a flow variable at the current slab are stored in the
F-array in a contiguous block or segment, as follows:

Two kinds of information are therefore needed in order to access these values:

- The starting position of the segment within the F-array
- The order in which the slab values are stored within the segment

The Zero-Location Index

The location just before the segment is called the zero-location index.

These zero-location indexes are allocated by PHOENICS at the beginning of the computation. Their values depend on the problem dimensionally, the variables solved for and stored, whether the problem is steady or transient, whether it is parabolic or elliptic, etc.

The integer function L0F (el-zero-ef) provides, given a flow-variable index, the zero-location index for the current slab of values.

For instance, the zero-location index for the values of W1 on the current slab is L0F(W1). The slab values of W1 are therefore stored in the segment:

F(L0F(W1)+1) to F(L0F(W1)+NX*NY)

The Segment Structure

The slab values for each variable are written in each segment in the following
order:

First, the NY cells at IX=1, then, the NY cells at IX=2 ..... finally, the NY cells at IX=NX.

Therefore, the cell (IX, IY) of the current slab of W1 is stored as

F(L0F(W1)+IY+NY*(IX-1)

The coding for the example in this lecture using direct access to the F-array would be:

196 CONTINUE C * ----------------SECTION 6----------------FINISH OF IZ SLAB L0W1 = L0F(W1) L0C1 = L0F(C1) L0C2 = L0F(C2) [ or L0C2 =L0F(LBNAME('MASS'))] DO IX=1, NX DO IY=1, NY I=IY+NY*(IX-1) F(L0C2+I)=F(L0W1+I)*F(L0C1+I) END DO END DO RETURN

Note that the zero-locations are obtained outside the loop, as they are invariant for the slab. They must be obtained at each slab, however, as they change from slab to slab.

The L0F function should never be used within an array reference, as this is very inefficient.

Note also that the double loop (IX, IY) above is in fact equivalent to the following single one:

DO I=1,NX*NY F(L0C2+I)=F(L0W1+I)*F(L0C1+I) END DO

As long as the loop covers the whole slab, this is more efficient than a double loop.

The Two Methods: Advantages And Disadvantages

Choosing one of the two methods is, to a great extent, a question of personal preference. However, these are some [subjective] advantages and disadvantages:

- m1: The FN library is a compact way of coding, which circumvents the use of do-loops. Furthermore, the FN subroutines have built-in checks that decrease the chances of an error. However, the number of FN functions available is, although large, limited, and the user might not find there what they are looking for. Computational efficiency is good.
- m2: The direct access to the F-array is a powerful, flexible and efficient method. However, there is no "safety net" once the zero-locations have been obtained; and changing the value of the wrong position in the F-array is a (very) difficult error to trace.

Users will probably find themselves moving from the first option to the last as they build up expertise.

The integer indices GRSP1 to GRSP10 provide access to auxiliary slabwise working space, which can be used in the calls to FN subroutines and in the F-array method.

Before one of these indices can be used, space must be allocated for it in the F-array.

This is effected in PHOENICS by a call to the subroutine MAKE in Group 1 Section 1.

Example:

C--- GROUP 1. RUN TITLE AND OTHER PRELIMINARIES C 1 GO TO (1001, 1002), ISC 1001 CONTINUE CALL MAKE(GRSP1) CALL MAKE(GRSP7)

Example of usage: Average value, at each cell, of the flow variables C1 to C10.

CALL MAKE(GRSP2) ! This in Group 1 Section 1 C... FN1(Y,A)..........Y=A ! This in, eg, Group 19 CALL FN1(GRSP2,0.0) ! Initialise GRSP2 to 0.0 DO IVARIA=C1, C10 C... FN60(Y,X)..........Y=Y+X CALL FN60(GRSP2, IVARIA) ! GRSP 2 will hold the sum of C's END DO C... FN25(Y,A)............Y=A*Y CALL FN25(GRSP2, 1.0/10.0) ! Now average

A similar device is frequently used in GREX, where another suite of indices (EASP1 to EASP20) are available. The EASP set can also be used in GROUND, but then GREX should be checked for conflicting use.

Fluid properties can be set in Group 9 of GROUND following the methods outlined above for flow variables.

The only difference is, in fact, the integer index used to refer to the property.

Two related quantities are used in GROUND in connection with properties. These are:

- The value of the property (e.g., RHO1), which is either a constant (RHO1=1.0) or a flag (if the property is to be computed in GROUND, e.g., RHO1=GRND).
- The integer index for a property, which is available for use even when the property is not STOREd in the Q1 file. The integer index is stated for each property at the top of its Section in Groups 9 and 10 (e.g., DEN1 for RHO1).

Storing properties is required for their relaxation and plotting with PHOTON. The SATELLITE STORE command, e.g., STORE(RHO1) or STORE(DEN1). EARTH will only allow access to the F-array via the correct integer indices.

A list of properties and indices can be found by looking at group 9 of the GROUND subroutine.

The properties - the thing set in Q1 to gain access to GROUND - are given on the left, the indices - used to set the properties into the F-array - are given on the right.

Example Of Setting A Property

We want to compute the temperature of the fluid as a function of enthalpy,
given by T=H1/Cp, where Cp is a constant. (Note that this option is already provided in
PHOENICS.)

Q1 settings:

TMP1=GRND; RG(1)=1.0E3

This instructs PHOENICS to compute the temperature in GROUND, and sends 1.0E3 as RG(1) to be used as Cp.

As a property is to be set, there is no question as to where the coding should be placed - it MUST be in Group 9, in the relevant section.

GROUND coding using the FN library:

900 CONTINUE C * ------------------SECTION 10------------------- C For TMP1.LE.GRND phase 1 temperature Index TEMP1 IF(TMP1.EQ.GRND) THEN C..........FN2(Y,X,A,B)..........Y=A+BX CALL FN2(TEMP1, H1, 0.0, 1.0/RG(1)) ENDIF RETURN

Note the IF....ENDIF protection. The section will be visited if TMP1 is set to any GRNDn number. The protection will prevent the unwanted execution of the coding in GROUND if TMP1 is not set to GRND.

Setting Boundary Conditions And Sources

Boundary conditions and sources are set in GROUND in Group 13.

When a GRNDn flag is used in a COVAL command as coefficient or value, then a particular section of Group 13 is visited by EARTH. EARTH expects the corresponding coefficient or value to be set in that section.

Coding can be inserted by the user following the same practices described for the flow variables. The integer index to be used is CO, if the coefficient is to be set, and VAL, if the value is to be set.

Variables are retrieved and set slabwise, as for the flow variables. However, when the FN library is used, PHOENICS automatically limits the changes introduced in the algebraic equations to those cells within the slab which are in the patch.

The name of the patch currently being dealt with is available in GROUND (character*8-variable NPATCH). This can be used in an IF...ENDIF structure as a further means of controlling the access to the coding.

This is the practice adopted in GREX, and its adoption is highly recommended.

The integer index of the flow variable being dealt with is also available as INDVAR, and can be used for flow control.

Example Of Source Setting

A chemical reaction is taking place in the computational domain, in which C1
and C2 produce C3:

C1 + C2 ---> C3.

The volumetric rate at which C1 is consumed is given by

rate = -A * C1 * C2 * T

where A is a constant and T is the temperature.

How can the sink of C1 be set up in PHOENICS? Three steps are needed:

- Expressing the source/sink in a linear form
- Providing the appropriate Q1 settings
- Providing the appropriate GROUND coding

Linearizing The Source

Sources and boundary conditions have to be cast into a coefficient/value form
before they can be handled by PHOENICS.

S_{f} = C(V_{f} -
f_{p})

where C is the coefficient and V is the value.

In the example considered, the sink of C1 can be expressed as:

S_{C1} = A * C2 * T * (0.0 - C1)

where the coefficient is C = A * C2 * T and the value is V = 0.0

Q1 Settings

In the Q1 file, PATCH and COVAL commands must be provided for this source, and
the value of A has to be transmitted to GROUND.

PATCH command:

PATCH(SINKOFC1,VOLUME,1,NX,1,NY,1,NZ,1,LSTEP)

Note that the patch covers the whole of the domain, and that the sink will be expressed per unit volume. (We are assuming that the reaction rate is a volumetric one).

COVAL command:

COVAL(SINKOFC1, C1, GRND, 0.0)

Since the coefficient is A * C2 * T, it has to be computed in GROUND. The value is, however, constant (V=0.0) and can be set in the COVAL command.

Transmission of A: The setting RG(1)=desired-value will suffice.

Ground Coding

__The coefficient C=A*C2*T has to be coded in GROUND. Since GRND is the flag used
in the coefficient slot of the COVAL command, the section of Group 13 to be visited is
Section 1. EARTH expects the coefficient for the current IZ slab to be set there: __

130 CONTINUE C-----------------SECTION 1---------------- coefficient = GRND IF((NPATCH.EQ.'SINKOFC1').AND.(INDVAR.EQ.C1)) THEN C.........FN21(Y, X1, X2, A, B) .....Y=A+B*X1*X2 CALL FN21(CO,C2,TEMP1,0.0,RG(1)) ENDIF RETURN

Note in the coding above:

- The use of the PATCH name and INDVAR for protection
- The use of CO as integer index for the coefficients

Note On Setting Boundary Conditions And Sources

As noted above, EARTH expects the user to set, in Group 13, the coefficients or
values for the current variable INDVAR at the current slab IZ.

But, of course, the PATCH considered might cover only a portion of the current slab, and not the whole of it.

The variables IXF, IXL, IYF, IYL (available in GROUND) are the first and last indexes of the cells in the current slab which are also in the current PATCH.

When CALLed from Group 13, the FN subroutines act over PATCHes (i.e., from IXF to IXL and from IYF to IYL), and not over the whole slab.

Z-direction control is provided automatically by EARTH - if the PATCH does not exist on
the current slab, GROUND will not be called for that PATCH.

Access To Geometrical Quantities

Integer indices are also provided for the access to cell areas, cell volumes and porosity factors.

These indices are:

- AEASTfor the cell East-area
- ANORTHfor the cell North-area
- AHIGHfor the cell high-area
- VOLfor the cell volume
- EPORfor the East-face porosity
- NPORfor the North-face porosity
- HPORfor the high-face porosity
- VPORfor the volume porosity

These indices can be used just as the flow-variable indices or the property indices discussed above.

In-Slab Distances

The following indices, which refer to x-direction distances for the current
slab, are also available in GROUND:

- XG2D addresses the x-coordinates of the cell centers;
- XU2D addresses the x-coordinates of the cell faces;
- DXG2D addresses the distances between adjacent cell-centers;
- DXU2D addresses the distances between adjacent cell-faces;

The corresponding indices are also available for the y-direction: YG2D, YV2D, DYG2D, DYV2D. For polar grids, two additional distances are available in the y (i.e., radial) direction:

RG2D radial distance to the cell center (taking RINNER into account)

RV2D radial distance to the cell face (taking RINNER into account)

All the indices in the previous slide must be allocated space in the F-array if they are to be used.

This is done through a call to the subroutine MAKE in Group 1 Section 1 for the indices
involved:

E.g.: CALL.MAKE(XU2D)

After this CALL, the indices can be employed in the usual way.

Out-Of-Slab Distances

__Z-wise distances can be also accessed in PHOENICS, through the following
integer indices:
__

- ZGNZ for the z-coordinates of the cell centers;
- ZWNZ for the z-coordinate of the cell faces;
- DZGNZ for the distance between cell centers;
- DZWNZ for the distance between cell faces.

These quantities can only be accessed through the subroutine GETZ (quantity index, local array, dimension of array) or through the L0F function.

Note that, for BFC's the access to distances, areas and volumes is through a special
subroutine, GTIZYX. See the Encyclopedia for a full description. These quantities may be
referenced directly, using L0B(variable) to obtain the zero-location indices.

Service Functions And Subroutines

PHOENICS provides a number of "service functions" and "subroutines" to promote economy and robustness in the coding.

The FN library , L0F and LBNAME are just three examples of them; some more will be given now.

Full lists of available subroutines and functions are given in the Encyclopedia, under Function subroutines.

Inspection of GREX and the GX routines will show many of these being used extensively.

Access To Neighboring Slabs

Three integer functions are provided for the access to neighboring slabs:

- HIGH(index) for the slab next to the current one
- LOW(index) for the slab before the current one
- OLD(index) for the same slab at the previous time-step

These integer functions take, as argument, an integer index (e.g., W1, or AEAST, or EASP2, or DEN1); and return another integer index which can be used for access to the corresponding slab of values.

For example, the following piece of coding would fill the slab values of C1 with the
previous-slab values

C..... FN0(Y,X).....Y=X CALL FN0(C1,LOW(C1))

Print-Out And Dump Subroutines

A set of subroutines (of generic name WRIT**) are provided for labeling the
print-out of integers, reals, logicals and characters.

The subroutine PRN will print the current-slab values of any EARTH variable which is
accessed through an integer index.

e.g. CALL PRN('w1 ',W1)

The subroutine DUMPS will produce, when called from GROUND, a file of all the EARTH arrays, in PHI-file format.

Flow-Control Functions

A number of logical functions are also provided, which can be used in IF...ENDIF blocks for flow control. Some examples are:

- The logical function FLUID1(variable-index), which is true if the variable is a first-fluid one
- The logical function STORE(variable-index), which is true if the variable is stored
- The logical function SOLVE(variable-index), which is true if the variable is solved for
- The logical function WHF(variable-index), which is true if the variable is to have its equations solved whole field

These are the main points conveyed in this lecture:

- Most of the PIL variables which are used in the Q1 file are also available in GROUND. They are, automatically, given the same values set by the user in the Q1 file (or their defaults), and can be used in any group of GROUND.
- PHOENICS provides also access to cell values of:
- flow variables;
- fluid properties;
- coefficients and values in source and boundary conditions;
- geometrical quantities.

- The key to access these quantities is always an integer index, which is used by PHOENICS to locate the variable within the F-array.
- The user can access the values of these quantities by two different methods:
M1: The FN library of subroutines, which perform operations with slab values directly on the F-array, in a manner which is fully transparent to the GROUND user

M2: Direct access to the F-array, in which the user modifies directly the values in the F-array

- PHOENICS provides also a tool-kit of service subroutines, intended to reduce the amount of coding users need to write, and therefore to reduce the possibility of errors.