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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.