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:
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 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:
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:
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:
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
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)
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 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:
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:
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:
Linearizing The Source
Sources and boundary conditions have to be cast into a coefficient/value form before they can be handled by PHOENICS.
Sf = C(Vf - fp)
where C is the coefficient and V is the value.
In the example considered, the sink of C1 can be expressed as:
SC1 = 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:
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:
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:
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:
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:
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:
These are the main points conveyed in this lecture:
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