2. Initial Project Study

As per the revised work plan, the first phase is of four months duration with a start date of January 1, 2009 and end date of April 30, 2009. This also included an interim report with a description of this phase along with planned code modifications and tests to be carried out in the next phase. There was no implementation work planned for this phase as it was to gain CITCOM related knowledge which was necessary in order to be able to work with the source code and then perform modifications and further implementations in the following two phases.

Four months duration of the first phase for the ``Initial Project Study'' is necessary for reasons explained in the subsequent paragraphs.

2.1 Background

This CITCOM package is PI's own copy of the Cartesian version of CITCOM source code and originated from the time he was working as postdoctoral researcher at the University of Colorado at Boulder, USA within Department of Physics during 2002-2004. In contrast to the Spherical version of the CITCOM which is well maintained and documented and is available at the Computational Infrastructure for Geodynamics web site [5], this Cartesian version has less documentation for the new users/developers to get kick started. Also, the source code is commented relatively sparsely otherwise our effort to get documentation extracted from the source code using Doxygen [8] would have been reasonably rewarding.

2.2 Learning Curve

Despite some difficulties related to understanding the source code due to relatively less documentation and sparsity of comments within the source code, getting to grips with the code was achieved by a number of ways including, but not limited to:

• Browsing through the source code to read the source code itself and sparse comments.
• Making use of Doxygen to extract as much as we could, at least "Call" graphs and "Called by" graphs were of help.
• Making use of etrace package [6] to learn function call tree.
• Meetings and discussions with CITCOM project PI, where his little-yellow-magic-notebook proved a very valuable asset.
• Googling.

As a result of our combined efforts, we became close to the point where with some confidence we could say that we are ready for the next phase of the project and to undertake code development work.

Efforts to learn more about the CITCOM package and implementation continue throughout this project. Getting hands on with the CITCOM packages, modifying existing source code, adding new functions, etc. is all useful.

2.2.1 Learning CITCOM

This CITCOM package is built on a structured finite element mesh which is made up of rectangular elements in two dimensions (2D) and brick elements in three dimensions (3D). These elements reduce to a square in 2D and a cube in 3D respectively. Given that an appropriate number of elements are used in each dimension with respect to the length of the corresponding dimension and as the size of each element along any axis is dependent of the number of elements and the overall length of the mesh along that direction. Such elements can be represented by figure 2.1 in the two dimensional case and by figure 2.2 in the three dimensional case respectively.

Figure 2.1: A standard rectangular element in 2D

Figure 2.2: A standard brick element in 3D

Our main focus here are three dimensional brick elements which for our application are more appropriate for the study of the convection of the earth's mantle. However, to make the learning systematic and smooth, we would continue with the description in two dimensions as well as in three dimensions. It is interesting to learn that structured mesh elements used by CITCOM, in contrast to the standard Cartesian coordinate system (x,y) in 2D and (x,y,z) in 3D, have a coordinates system (x,z) in 2D and (x,z,y) in 3D as shown in figure 2.3 for the two dimensional rectangular element and in figure 2.4 for the three dimensional brick element respectively.

Figure 2.3: A CITCOM rectangular element in 2D

Figure 2.4: A CITCOM brick element in 3D

This fact is further highlighted in figure 2.5 for the two dimensional rectangular element and in figure 2.6 for the three dimensional brick element along with given coordinates for a unit element where it is also reflected that the values of z-coordinates along the downward z-axis are taken positive.

Figure 2.5: A CITCOM rectangular element and local nodes relationship

Figure 2.6: A CITCOM brick element and local nodes relationship

Local node numbering for each element is counter-clockwise starting at the origin (0,0) in 2D as shown in figure 2.5. This is further elaborated in the 3D case where the local node numbering for an element starts at the origin (0,0,0), first on the front-face of the element followed by on the back-face of the element as shown in figure 2.6. This local node numbering runs from 1 to 4 in 2D and from 1 to 8 in 3D with the 0 (zero) position leaving un-used, in contrast to the C language standard (the language of the CITCOM package) where everything begins at 0. This is a major feature of CITCOM that almost everywhere the 0 (zero) position is un-used and an extra position is allocated where needed in lieu of this un-used position. However, there are a few exceptions.

There is another point which could be of interest. In most cases, instead of one extra position in lieu of un-used 0 (zero) position, two extra positions are defined/allocated.

Also, no element on the mesh surface is considered as a surface element except those which are on the earth's surface considering the mapping between mesh and earth. In other words, only those elements which touch the X-axis in two dimensions and XOY-plane in three dimensions are considered to be surface elements. Similarly, only nodes which are on the X-axis in two dimensions and XOY-plane in three dimensions are considered to be surface nodes.

To elaborate further on this, in the two dimensional case as shown in figure 2.7, element 4 is not a surface element whereas element 7 is a surface element. Also, element 4 has no surface node whereas element 7 has two surface nodes which are 10 & 13. To emphasise this, any element with surface nodes will have exactly two nodes on the surface irrespective of element position within the mesh.

Figure 2.7: A sample CITCOM mesh of 8 (4x2) rectangular elements

Similarly, in the three dimensional case as shown in figure 2.8, element 4 is not a surface element whereas element 15 is a surface element. Also, element 4 has no surface node whereas element 15 has four surface nodes which are 40, 25, 28 & 43. To emphasise this, any element with surface nodes have exactly four nodes on the surface irrespective of element position in the mesh.

Figure 2.8: A sample CITCOM mesh of 16 (4x2x2) brick elements

It should be noted that the previous paragraphs which describe surface elements in a mesh are from the "numerical" surface of a finite element mesh perspective rather then from the "earth science" perspective. In the later case, the scenario is obvious and well understood. Hence this might just be a nomenclature issue but it is included here for the sake of clarity.

Mapping of local node numbering for any given element to the corresponding global node numbering is another important relationship for any finite element mesh. This helps to establish mesh connectivity and the relationship between the nodal coordinates for any given node in the mesh. Such a connectivity relationship is depicted in figure 2.7 for element 4 and element 7 for the two dimensional case and in figure 2.8 for element 4 and element 15 for the three dimensional case. Having established these relationships for a finite element mesh, it is not hard to compute the required parameters of the mesh such as global nodes for other given parameters such as elements and vice versa. For example, for any given element number, it is not hard to find a global node number corresponding to a local node number for that element and hence coordinates for the global node. Conversely, for any given global node number it can easily be established just how many elements in the mesh this particular node belongs to. In all mesh elements, local nodes, global nodes, node coordinates, relationship, etc. extra complexity is introduced by the levels of mesh and all of these entities exist at each level of the mesh. Each level in the mesh represents a certain depth of when an element is refined or sub-divided, for example, if a certain element in the mesh is refined or sub-divided three times, all the resulting elements belong to the third level of the mesh. These levels of the finite element mesh play an important role in the multigrid algorithms implemented (or to be implemented) in the CITCOM package. Rather, it is more appropriate to say that these levels in the finite element mesh provides the basis for multigrid algorithms.

Let us come back to the mesh element-node relationship. We notice that all numbering for elements, local nodes and global nodes starts at 1 rather than 0 (zero). This element numbering and global node numbering begins at the origin of the mesh, and gets incremented along the z-axis, followed by the x-axis and in the 3D case followed by the y-axis. This incremental order for the element and global node numbering plays an important role in calculating any entity, whether it is an element-node relationship or the computation of velocity, pressure, viscosity, etc.

In addition to this order of element and global node numbering in the mesh, a number of C data structures defined in the header files are of great importance in implementing algorithms in CITCOM but probably it is not sensible to go into those details here.

2.3 Precautions

There is no restriction on the number of processes any CITCOM job can run for or any specific sequence of the number of processes such as 2, 4, 8, .... However, there are some very basic rules which must be met for any CITCOM job to be completed successfully. In this regard, a few precautions, as described below, must be taken.

• The number of mgunit in the (x, z, y)-direction must be at least twice the number of processes in the (x, z, y)-direction respectively. These mgunit<x,z,y> define the number of elements along the (x, z, y)-direction for the corresponding structured finite element mesh at the coarsest level. Thus, the exact number of elements and the number of global nodes in the (x, z, y)-direction for the given mesh level are defined as:
  • number of elements = mgunit<x, z, y> * $2^{(level-1)}$
  • number of nodes = mgunit<x, z, y> * $2^{(level-1)} + 1$

  This helps determine the number of global nodes along the (x, z, y)-directions to be specified in the CITCOM input file.

  The description for the three dimensional case as given above is also valid for the two dimensional case.

• Any number of processes can be defined in the CITCOM input file as long as they are consistent with the other predefined parameters, except the number of processes along the z-direction which must be 1 (one).

2.4 Miscellaneous

This section includes the work which may not necessarily contribute towards the improvement of the performance of CITCOM but might improve the CITCOM code in order to make it more user friendly and easy to understand. The first of these such additions is described below, others could be added later as and when implemented.

2.4.1 data_directory()

The CITCOM code reads in all of its required parameters from an input file which is read by CITCOM as first and only input argument at the command-line. In this input file, a user is required to provide an output directory and the base name for the output files in the form of <dirname>/<basefilename> against the input keyword datafile. (Note: The user is free to add new keywords and/or change existing keywords. Precautions must be taken as this may trigger changes to the code itself unnecessarily.)

In order to run CITCOM, a user is required to create an empty output directory before submitting the job, or more appropriately before the start of execution, within the directory from where the CITCOM job will run. If the user forgets to create the directory and the job starts executing, it is certain that the job will crash. However, if the directory was specified in the input file as existing from a previous run, every thing will be overwritten unless the base file name is different for the current run.

Addition of the function data_directory() has eliminated this unpleasant possibility. The iser is no longer required to create a directory as specified in the input file. It will be created at run time. If the directory with the same name as specified in the input file exists, it will be renamed to a directory, the name of which is generated by attaching the current date and time (as explained in section 2.4.2) with current process identity as an extention to the existing name in the form <directory>.<yyyymmdd>-<hhmmss>.<pid> to make it unique. Then, the directory as specified in the input file is created by the root or 0 process and it is assured that no process is allowed to write to any of the output files until this new directory is created.

This would ensure that no job ends up wasting compute resources and no results from a previous run are over written due to minor negligence at the user end.

2.4.2 date_and_time()

This function utilises C' built in time struct (struct tm) and returns the current date and time in the format <yyyymmdd>-<hhmmss> as required by the function date_directory() (see section 2.4.1) to make the output data directory name unique.

2.5 Planned Implementation and Tests

2.5.1 Implementation

During the second phase of this dCSE project, the implementation tasks as defined in the project proposal are to be implemented. These tasks are outlined very briefly as follows:

• Study of the currently implemented multigrid V-cycle.
• Implementation of multigrid W-cycle.
• Implementation of multigrid F-cycle.

2.5.2 Testing

Once the implementation of these proposed multigrid cycles is complete, we would go through testing all three multigrid cycle implementations against a comparatively smaller to medium size test problem. Hence the tasks are:

• Test multigrid V-cycle for a smaller to medium size test problem(s).
• Test multigrid W-cycle for a smaller to medium size test problem(s).
• Test multigrid F-cycle for a smaller to medium size test problem(s).

2.5.3 Optimisation

All the test results will then be compared against standard benchmarks for the selected test problem(s). After verifying the accuracy of the test results, the implementation of multigrid cycles will go through an optimisation phase. This will require completion of the following tasks:

• Optimise multigrid V-cycle implementation.
• Optimise multigrid W-cycle implementation.
• Optimise multigrid F-cycle implementation.

At this point it would be especially interesting to see/note if the optimisation would be problem-dependent, and if so, what characteristics determine the optimal conditions for these cycles.

2.5.4 Further Testing

Next comes the testing of optimised multigrid cycles for test problem(s) ranging from medium to large and possibly very large size and subject to the feasibility of any such test problem in the context of CITCOM capability and the limitations of the problem itself.

• Test multigrid V-cycle for medium to large size test problem(s).
• Test multigrid W-cycle for medium to large size test problem(s).
• Test multigrid F-cycle for medium to large size test problem(s).

After satisfactory completion of the test runs and verification of results, other aspects of the CITCOM package would be considered for optimisation as appropriate. This should lead to the next and final test runs for a number of benchmark test problems by the end of which we are expecting to have reasonable performance results available for the report on the end of the second phase.

2.6 Conclusions

On completion of the first phase of the project within prescribed timescale, multigrid cycles implementation will be undertaken in the next phase. On the basis of lessons learnt, it has been understood and recognised mutually that there should have been more time dedicated towards the learning and better understanding of this package. In line with this need, time will be spared to learn more about the CITCOM source code during the next two phases on the basis of as and when feasible. At the same time, it been agreed that we will continue to document and improve the existing documentation both in the form of source code commentary and this report.