## 4 Theory and Calculation
### Temporary Array Method  
A temporary array-based technique can be used to ensure reproducibility for incrementing operations. Consider using an edge-\(\rightarrow\)cells mapping and an incrementing operation. Here, we would iterate through all the edges, calculate values, and add them to a variable defined on a neighboring cell. To achieve reproducibility, we modify this structure by storing the calculated increments in a temporary array defined on the edges, and after all the increments are calculated, we iterate through all the cells and apply these increments in a fixed order defined by the global_IDs of the edges. In Figure 5 we can see an example of this method, where edge2's global_ID is the smallest, so the value from edge2 is applied first on the cell, then edge0, etc. To achieve this modified execution, a few extra preparations must be conducted in the backend, which are shown in Algorithm 1. After the global_IDs are shared, the next step is to create a reversed mapping for every map. The reversed mapping is needed so we can iterate through the cells and in each iteration we can access the edges connected to the given cell. This reversed map uses local indices which might be in different order in different MPI ranks. That is why we need to reorder them by using their previously shared global indices. Another modification conducted on the reversed map is that it actually stores indices of a temporary array where the increments from the edges are stored for a cell. In other words, if the \(k\)th element in line \(n\) (\(k\)th edge connection of cell \(n\)) of the  
Figure 5: Example execution order of edges around a cell. Due to local id renumbering, the global ids must be used for a reproducible execution order. reversed map is \(x\), then it means that in the temporary increments array at location \(x\) the increment for cell \(n\) from edge \(k\) can be found. The main disadvantage of this method is the need for significant additional memory to store the reversed mapping, and to store the increments. The reversed map uses a Compressed sparse row (CSR) format, which consists of a main array of increment indexes (integers), with the size of set_from_size * original_map_dimension, and another array indexing the previous array with a size of set_to_size\(+1\). The temporary arrays themselves can use much more memory: set_from_size * map_dimension * data_element_size. ```
exchange global IDs \(OP\_map\_index\) = number of maps for\(m=0\) to \(OP\_map\_index\)do  create reversed mapping for map \(m\) \(set\_to\_size\) = target set's size of map \(m\) for\(i=0\) to set_to_size do  sort the reversed connections of \(i\) by global IDs endfor endfor
```  
**Algorithm 1** Algorithm of generating incrementing order  
After creating the reversed map with the correct order, we generate a new op_par_loop implementation code to use this modified method. The main changes can be seen in Algorithm 2. After the initialization phase, it is imperative to set all elements in a temporary array to zero to accommodate individual increments. This step is crucial as the user kernel performs the increments, and proper initialization is required beforehand. Moreover, this approach ensures that the data remain in the cache, enhancing overall performance. Then we can call the kernel function for all edges to access the elements defined on the cells. If a parameter is accessed through an OP_READ or OP_WRITE method, then the execution order does not matter, so we can use the original method of directly storing the new state in the data. If the parameter is incremented (OP_INC), then we need to store each increment value in the tmp_incs array instead of adding it to the actual data. After the iteration on edges is completed and all increments are calculated, we need to apply those to the actual data on cells. For that, we start a new cell-based loop on the cells and by using the reversed mapping with the fixed ordering, for each cell, we can gather and apply the increments.