Replacing For-Loops with a Macro in C functions

[mabruzzo]

This issue documents a way to potentially enhance the codebase that briefly discussed in PR #106.

Overview

The idea is to replace nested for-loop in (C/C++) code with a macro of the format: FULL_LOOP(index, index_helper, omp_pragma), where

• index is the name of the variable that is being iterated over
• index_helper is a previously initialized instance of grackle_index_helper
• omp_pragma is a string-literal holding the OpenMP-related pragma arguments. This is conditionally applied (when _OPENMP is defined) to the outer-loop using the C99 _Pragma keyword.

A major motivation for doing something like this is to facilitate the testing of different loop structures for GPU parallelism in the future.

For concreteness, I'll show how one of the loops in local_calculate_pressure would change.

Effective current implementation (with changes from PR #106)

  double tiny_number = 1.e-20;
  const grackle_index_helper ind_helper = _build_index_helper(my_fields);

# ifdef _OPENMP
# pragma omp parallel for schedule( runtime )
# endif
  for (int outer_ind = 0; outer_ind < ind_helper.outer_ind_size; outer_ind++){

    const grackle_index_range range = _inner_range(outer_ind, &ind_helper);

    for (int index = range.start; index <= range.end; index++) {

      pressure[index] = ((my_chemistry->Gamma - 1.0) *
			 my_fields->density[index] *
			 my_fields->internal_energy[index]);
 
      if (pressure[index] < tiny_number)
        pressure[index] = tiny_number;
    } // end: loop over i
  } // end: loop over outer_ind

Equivalent implementation with FULL_LOOP macro

  double tiny_number = 1.e-20;
  const grackle_index_helper ind_helper = _build_index_helper(my_fields);

  FULL_LOOP(index, ind_helper, "omp parallel for schedule( runtime )") {

    pressure[index] = ((my_chemistry->Gamma - 1.0) *
                       my_fields->density[index] *
                       my_fields->internal_energy[index]);
 
    if (pressure[index] < tiny_number)
      pressure[index] = tiny_number;
  }

Sample Implementation

An example implementation of this macro is provided below (I have confirmed that implementation works - with & without OpenMP):

#ifdef _OPENMP
#define APPLY_OMP_PRAGMA(statement) _Pragma(statement)
#else
#define APPLY_OMP_PRAGMA(statement) /* ... */
#endif

#define GET_J_(outer_index, index_helper)                             \
  ((outer_index % index_helper.num_j_inds) + index_helper.j_start)
#define GET_K_(outer_index, index_helper)                             \
  ((outer_index / index_helper.num_j_inds) + index_helper.k_start)

#define FULL_LOOP(index, index_helper, omp_pragma)                        \
  APPLY_OMP_PRAGMA(omp_pragma)                                            \
  for (int oind_ = 0; oind_ < ind_helper.outer_ind_size; oind_++)         \
    for (int iind_ = index_helper.i_start,                                \
         index = (index_helper.i_start + index_helper.i_dim *             \
                      (GET_J_(oind_, index_helper) + index_helper.j_dim * \
                       GET_K_(oind_, index_helper)));                     \
         iind_ <= index_helper.i_end; iind_++, index++)

I'm not a huge fan of declaring/incrementing 2 variables inside the inner-loop, but doing it this way makes the application of the loop a lot cleaner. I have confirmed that the above macro definitely works (both with and without openmp), by converting both loops in local_calculate_pressure.

At this point, it would not take much effort to finish implementing this change, so let me know if you want me to do that (in PR #106 or a separate PR). However, it may make more sense to wait until we really need this for GPU-experimentation