More OpenMP usage in equations and integrators.

Earlier we only used OpenMP for the Equation's `loop` and `loop_all`
methods.  This is not optimal as the number of cores increase.  With
12-20 core machines being routinely available it makes sense to use
OpenMP as much as possible and so now all the equation methods,
`initialize`, `initialize_pair`, `post_loop` in addition to the original
`loop/loop_all` are executed in parallel.  In addition, the integrator
steps are also executed in parallel.  This should give us much better
scale up with more cores.  Note that in the code, the loops do not set
nogil=True since they are already in a block where the gil is released.
This PR requires the latest version of compyle with
pypr/compyle#37 merged.

pysph/sph/acceleration_eval_cython.mako

@@ -37,7 +37,7 @@ ${indent(all_eqs.get_py_initialize_code(), 0)}
 #######################################################################
 % if all_eqs.has_initialize():
 # Initialization for destination ${dest}.
-for d_idx in range(NP_DEST):
+for d_idx in ${helper.get_parallel_range("NP_DEST")}:
     ${indent(all_eqs.get_initialize_code(helper.object.kernel), 1)}
 % endif
 #######################################################################
@@ -46,7 +46,7 @@ for d_idx in range(NP_DEST):
 % if len(eqs_with_no_source.equations) > 0:
 % if eqs_with_no_source.has_loop():
 # SPH Equations with no sources.
-for d_idx in range(NP_DEST):
+for d_idx in ${helper.get_parallel_range("NP_DEST")}:
     ${indent(eqs_with_no_source.get_loop_code(helper.object.kernel), 1)}
 % endif
 % endif
@@ -65,7 +65,7 @@ ${indent(helper.get_src_array_setup(source, eq_group), 0)}
 src_array_index = src.index
 
 % if eq_group.has_initialize_pair():
-for d_idx in range(NP_DEST):
+for d_idx in ${helper.get_parallel_range("NP_DEST")}:
     ${indent(eq_group.get_initialize_pair_code(helper.object.kernel), 1)}
 % endif
 
@@ -78,7 +78,7 @@ nnps.set_context(src_array_index, dst_array_index)
 ${helper.get_parallel_block()}
     thread_id = threadid()
     ${indent(eq_group.get_variable_array_setup(), 1)}
-    for d_idx in ${helper.get_parallel_range("NP_DEST")}:
+    for d_idx in ${helper.get_parallel_range("NP_DEST", nogil=False)}:
         ###############################################################
         ## Find and iterate over neighbors.
         ###############################################################
@@ -105,7 +105,7 @@ ${helper.get_parallel_block()}
 ###################################################################
 % if all_eqs.has_post_loop():
 # Post loop for destination ${dest}.
-for d_idx in range(NP_DEST):
+for d_idx in ${helper.get_parallel_range("NP_DEST")}:
     ${indent(all_eqs.get_post_loop_code(helper.object.kernel), 1)}
 % endif
 

pysph/sph/acceleration_eval_cython_helper.py

@@ -272,8 +272,11 @@ def get_parallel_block(self):
         else:
             return "if True: # Placeholder used for OpenMP."
 
-    def get_parallel_range(self, start, stop=None, step=1):
-        return get_parallel_range(start, stop, step)
+    def get_parallel_range(self, start, stop=None, step=1, nogil=True):
+        if nogil:
+            return get_parallel_range(start, stop, step, nogil=True)
+        else:
+            return get_parallel_range(start, stop, step)
 
     def get_particle_array_names(self):
         parrays = [pa.name for pa in self.object.particle_arrays]

pysph/sph/integrator_cython.mako

@@ -102,7 +102,7 @@ cdef class Integrator:
         # Only iterate over real particles.
         NP_DEST = dst.size(real=True)
         ${indent(helper.get_array_setup(dest, method), 2)}
-        for d_idx in range(NP_DEST):
+        for d_idx in ${helper.get_parallel_range("NP_DEST")}:
             ${indent(helper.get_stepper_loop(dest, method), 3)}
         % endif
         % endfor

pysph/sph/integrator_cython_helper.py

@@ -15,6 +15,7 @@
 from pysph.sph.equation import get_array_names
 from .acceleration_eval_cython_helper import get_helper_code
 from compyle.api import CythonGenerator, get_func_definition
+from compyle.cython_generator import get_parallel_range
 
 
 getfullargspec = getattr(
@@ -66,6 +67,12 @@ def get_helper_code(self):
         code = get_helper_code(helpers)
         return '\n'.join(code)
 
+    def get_parallel_range(self, start, stop=None, step=1, nogil=True):
+        if nogil:
+            return get_parallel_range(start, stop, step, nogil=True)
+        else:
+            return get_parallel_range(start, stop, step)
+
     def get_stepper_code(self):
         classes = {}
         for dest, stepper in self.object.steppers.items():