Merge pull request #554 from teeteg/fast_secondaries

fix and enhancement to secondary propagation in ModuleList

Author: JulienDoerner

CHANGELOG.md

@@ -1,11 +1,13 @@
 ## CRPropa vNext
 
 ### Bug fixes:
+* Fixed the secondariesFirst parameter not actually allowing secondary particles to be propagated before the primary
 * threshold calculation in PhotoPionProduction
 
 ### New features:
 
 * Improved plugin-template and added test to ensure the template is working
+* Secondaries are now distributed over multiple threads if OpenMP parallelisation is enabled
 
 ### Interface changes:
 

include/crpropa/ModuleList.h

@@ -82,40 +82,48 @@ class ModuleList: public Module {
 	 * @param recursive  Whether to also process possible created secondaries
 	 * @param secondariesFirst  Whether to process the secondaries directly after the
 	 * step where they are created or to only process them after the primary is fully finished
+	 * @param waitForSecondaries  Whether to wait for the secondaries to be fully processed 
+	 * before continuing with the primary candidate if secondariesFirst is set to true
 	 */
-	void run(Candidate* candidate, bool recursive = true, bool secondariesFirst = true);
+	void run(Candidate* candidate, bool recursive = true, bool secondariesFirst = true, bool waitForSecondaries = true);
 	/** Run function
 	 * This function wrapps the ModuleList::run(Candidate* candidate) function
 	 * This run function does the full simulation for a single given candidate
 	 * @param candidate  Candidate to use for the simulation
 	 * @param recursive  Whether to also process possible created secondaries
 	 * @param secondariesFirst  Whether to process the secondaries directly after the
 	 * step where they are created or to only process them after the primary is fully finished
+	 * @param waitForSecondaries  Whether to wait for the secondaries to be fully processed 
+	 * before continuing with the primary candidate if secondariesFirst is set to true
 	 */
-	void run(ref_ptr<Candidate> candidate, bool recursive = true, bool secondariesFirst = true);
+	void run(ref_ptr<Candidate> candidate, bool recursive = true, bool secondariesFirst = true, bool waitForSecondaries = true);
 	/** Run function
 	 * This function calls the run function for a single Candidate on each entry of the given candidate_vector.
 	 * If OpenMP parallelization is activated (default) the execution of the run function for each member
-	 * of candidate_vector is parallelized. However, the secondaries are still processed in order.
+	 * of candidate_vector is parallelized. 
 	 * @param candidates  Candidate vector
 	 * @param recursive  Whether to also process possible created secondaries
 	 * @param secondariesFirst  Whether to process the secondaries directly after the
 	 * step where they are created or to only process them after the primary is fully finished
+	 * @param waitForSecondaries  Whether to wait for the secondaries to be fully processed 
+	 * before continuing with the primary candidate if secondariesFirst is set to true
 	 */
-	void run(const candidate_vector_t *candidates, bool recursive = true, bool secondariesFirst = true);
+	void run(const candidate_vector_t *candidates, bool recursive = true, bool secondariesFirst = true, bool waitForSecondaries = true);
 	/** Run function
 	 * This function generates count Candidates over the given source.
 	 * Each generated Candidate is then handed over to ModuleList::run(Candidate*)
 	 * where it is simulated.
 	 * If OpenMP parallelization is activated (default) the count primary Candidates are
-	 * generated and simulated in parallel. However, the secondaries are still processed in order.
+	 * generated and simulated in parallel. 
 	 * @param source  The source to generate Candidates from
 	 * @param count  The number of Candidates to simulate
 	 * @param recursive  Whether to also process possible created secondaries
 	 * @param secondariesFirst  Whether to process the secondaries directly after the
 	 * step where they are created or to only process them after the primary is fully finished
+	 * @param waitForSecondaries  Whether to wait for the secondaries to be fully processed 
+	 * before continuing with the primary candidate if secondariesFirst is set to true
 	 */
-	void run(SourceInterface* source, size_t count, bool recursive = true, bool secondariesFirst = true);
+	void run(SourceInterface* source, size_t count, bool recursive = true, bool secondariesFirst = true, bool waitForSecondaries = true);
 
 	/// @return Returns the description string for all modules
 	std::string getDescription() const;

src/ModuleList.cpp.in

@@ -61,40 +61,54 @@ void ModuleList::process(ref_ptr<Candidate> candidate) const {
 	process((Candidate*) candidate);
 }
 
-void ModuleList::run(Candidate* candidate, bool recursive, bool secondariesFirst) {
+void ModuleList::run(Candidate* candidate, bool recursive, bool secondariesFirst, bool waitForSecondaries) {
+	
+	int taskedSecondaries = 0; //Secondaries already submitted as tasks for all threads
+	
 	// propagate primary candidate until finished
 	while (candidate->isActive() && (g_cancel_signal_flag == 0)) {
 		process(candidate);
 
 		// propagate all secondaries before next step of primary
 		if (recursive and secondariesFirst) {
-			for (size_t i = 0; i < candidate->secondaries.size(); i++) {
+			#pragma omp taskloop num_tasks(2 * omp_get_max_threads()) untied
+			for (size_t i = taskedSecondaries; i < candidate->secondaries.size(); i++) {
 				if (g_cancel_signal_flag != 0)
-					break;
+					continue;
 				run(candidate->secondaries[i], recursive, secondariesFirst);
 			}
+
+			// increase counter of already tasked secondaries
+			taskedSecondaries += candidate->secondaries.size();
+			
+			// wait for completion of tasked secondaries or continue with primary
+			if (waitForSecondaries) {
+				#pragma omp taskwait
+				} 
 		}
 	}
 
 	// propagate secondaries after completing primary
 	if (recursive and not secondariesFirst) {
-		for (size_t i = 0; i < candidate->secondaries.size(); i++) {
+		#pragma omp taskloop num_tasks(2 * omp_get_max_threads()) untied
+		for (size_t i = taskedSecondaries; i < candidate->secondaries.size(); i++) {
 			if (g_cancel_signal_flag != 0)
-				break;
+				continue;
 			run(candidate->secondaries[i], recursive, secondariesFirst);
 		}
+		#pragma omp taskwait
 	}
 
 	// dump candidae and secondaries if interrupted.
 	if (candidate->isActive() && (g_cancel_signal_flag != 0)) 
 		dumpCandidate(candidate);
 }
 
-void ModuleList::run(ref_ptr<Candidate> candidate, bool recursive, bool secondariesFirst) {
-	run((Candidate*) candidate, recursive, secondariesFirst);
+void ModuleList::run(ref_ptr<Candidate> candidate, bool recursive, bool secondariesFirst, bool waitForSecondaries) {
+	run((Candidate*) candidate, recursive, secondariesFirst, waitForSecondaries);
 }
 
-void ModuleList::run(const candidate_vector_t *candidates, bool recursive, bool secondariesFirst) {
+void ModuleList::run(const candidate_vector_t *candidates, bool recursive, bool secondariesFirst, bool waitForSecondaries) {
 	size_t count = candidates->size();
 
 #if _OPENMP
@@ -122,7 +136,7 @@ void ModuleList::run(const candidate_vector_t *candidates, bool recursive, bool
 		}
 
 		try {
-			run(candidates->operator[](i), recursive);
+			run(candidates->operator[](i), recursive, secondariesFirst, waitForSecondaries);
 		} catch (std::exception &e) {
 			std::cerr << "Exception in crpropa::ModuleList::run: " << std::endl;
 			std::cerr << e.what() << std::endl;
@@ -150,12 +164,11 @@ void ModuleList::run(const candidate_vector_t *candidates, bool recursive, bool
 	}
 }
 
-void ModuleList::run(SourceInterface *source, size_t count, bool recursive, bool secondariesFirst) {
+void ModuleList::run(SourceInterface *source, size_t count, bool recursive, bool secondariesFirst, bool waitForSecondaries) {
 
 #if _OPENMP
 	std::cout << "crpropa::ModuleList: Number of Threads: " << omp_get_max_threads() << std::endl;
 #endif
-
 	ProgressBar progressbar(count);
 
 	if (showProgress) {
@@ -189,7 +202,7 @@ void ModuleList::run(SourceInterface *source, size_t count, bool recursive, bool
 
 		if (candidate.valid()) {
 			try {
-				run(candidate, recursive);
+				run(candidate, recursive, secondariesFirst, waitForSecondaries);
 			} catch (std::exception &e) {
 				std::cerr << "Exception in crpropa::ModuleList::run: " << std::endl;
 				std::cerr << e.what() << std::endl;