I am trying to set up a 2 Phase simulation scheme, since in our project we are using quite a thick shielding and i don’t want to use expensive simulation time on the interactions in the shielding every single time the geometry inside the shielding is changed. We are interested in the activation outside of the shielding though, so I need to take the outer fluences into account. For this purpose I have set up a User Routine employing MGDraws BDXDRAW to create the necessary Phase Space files for the second part of the Simulation. In order to get enough statistics I have used AutoBiasing in the past, but now we want to use Clusters to further increase the statistics. My Questions would be the following:
When using MGDraw to dump the information of the boundary crossing particles you can output the particle weight, but when this is used simultaneously with AutoBiasing, is there a problem i am currently overlooking? Because currently these weights are really low, at around 10^-12 to 10^-15. And furthermore can i simply combine the output of these files, when using severa spawns? Ofc taking into account that these have to be divided by the total number of primaries used. I have attached one of those output files to this post.
When using a large number of primaries, there are some error messages in the .err file, that hint to a mistake in the geometry tracking i.e. “emfGeo: Ncoun 2000000000” or “Emsnsc”. I am currently unsure what some of those errors mean, since I can’t seem to find an explanation of i.e. Emsnsc in the documentation or in the slides of the beginners course.
I have read in this forum that these kind of errors are not breaking the simulation, but ofc they are unwanted, because i would assume these errors are the reason for the increased time per primary:
What could be a good practice to debug the geometry regarding these issues. I know that it is recommended not to use touching bodies, but in a constantly changing geometry, due to the developement process, i found it rather complicated to only use plains, because the code gets unreadable rather fast.
When using a Cluster with a 24h limit to calculations, especially when using AutoBiasing, is it advantages to use one cylce and just a time limit without a primary limit? Currently I am using a several cycle structure with a limit in primaries and additionaly a time limit, so the run doesnt crash from the cluster side. This is due to the complications in calculating the expected time for a set amount of primaries due to the large discrepancy in CPU time per primary. But i guess, i only need the information in the .out file how large the total weight of primaries per run was, since the source is currently unbiased with primary weight 1.0.
I will go through your different questions point by point while trying to answer all the doubts.
Question 2a: low particle weight
I confirm that I can reproduce the same particle weights that you report. Looking at the heavy shielding and the spread ROI regions set in the AUTOIMBS card, it is not too surprising that the biasing scheme applied is quite strong and the resulting weights are low.
Question 2b: combination of multiple mgdraw phase space file
Yes, as you mention, you can combine the output files produced by different spawns while considering the total number of primaries (or better the total weight of the primary particles) simulated.
Question 3: error messages in the .err file and higher Maximum CPU time compared to Average *** Emfgeo: Ncoun, 2000000000 is a counter that is printed every time that (cumulatively) electrons, positron and photon do 2 billion steps, and so approaching the upper limit for an integer. This specific output-error is reported in the manual as for the Umfnst: one. *** Emsnsc: sighlp,sigmlr,eusplc is a printing done by the routine for single scattering of electrons and positrons. You can disregard this error, as with the others too.
The higher Maximum CPU time compared to the Average CPU time can be related to the use of AUTOIMBS. Given the learning-based nature of the tool, some particles are simulated with little or no biasing at the start, and therefore they can require a larger CPU time compared to the last ones.
Question 4: good practice to debug the geometry
Your simulation geometry looks error-less. Generally, in addition to the recommendation of using planes and not touching bodies, checking the error file to spot eventual geometry error is a good practice that you are already following.
Question 5: time limit for run on cluster
The number of primaries is the variable driving the uncertainty of the simulation. For this reason I think it is good to have direct control on it. Said so, you have the flexibility of running the simulation with a limit imposed by the time rather than a number of primary. The additional drawback I can think about is that you will always need to retrieve for each cycle run in the cluster the actual number of simulated particle (which is a quite scriptable task but could be annoying).
Just one additional note. In the mgdraw.f routine, you are writing out the value of ETRACK which contains the total energy of the particle (kinetic energy + mass).
thank you very much for the extensive answers! They clear up a lot of things.
It was always a bit confusing to look at the avg. and max. time per primary and have a difference of about 4 orders of magnitiude. I have written a small analysis programm for all the .out and .error files of all spawns on the cluster and it seems to be occuring in all of these runs, with the max. time ranging from 1700 sec to 5900 sec. That would also explain, why flair indicates in a local test run that the simulation is “timed out” when a single primary sometimes takes more then 1.5 hours.
But your answers seem to indicate that these long computational times mostly occure due to the AutoImbs and the extensiv shielding and are not a result of faulty geometrical definitions and numerical errors. Therefore i need other ways to improve the statistics, becaus eeven when combining the phase space files, the summed weight of all 4.6\cdot10^6 individual scored photons leaving the Housing is 3.36. Similarly the total weight of all 2.5\cdot10^6 neutrons is 4.48, whilst all runs combined have a total primary weight of 7.6\cdot10^6.
So the followup question would be (maybe that is a topic for another thread): Do I just throw more computational power and time at that problem?
Could another optimization strategy be, to employ a 3 Stage scheme without the AutoImbs function?
Simulating the inner region and scoring the particles entering the shielding.
Then using these phase space files to do the shielding calculations with a manual biasing scheme and scoring the leaving particles
Then using these phase space files for the outside calculations.
But this 3 phases probably aren’t that much faster and more prone to overlooking statistical mistakes.
The 3-steps approach you define can be useful to have at least a clear understanding of the energy spectra for the different particles that are going to be shielded (i.e the result of the first step). However, this might provide little optimization in terms of CPU-time since more than 50% of it is used to simulate the particles in the first shielding layer. This is shown from the first lines of the results given by the PROFILE (manual) card obtained running 1e4 primaries:
# Region Total(ms) % p n e+- g
0 -total- 60.34 100.00 1.5 0.8 17.5 78.8
3 L1Pb 32.48 53.82 0.0 0.5 0.7 98.4
19 TarFront 12.97 21.50 2.9 0.0 72.0 25.0
2 HouseIn 5.31 8.79 0.0 1.9 5.0 92.3
It is clear that the most time-consuming process is the transport of gammas. Therefore, increasing the EMF transport and PROD-CUT threshold (avoiding the description of a whole electromagnetic shower that might be irrelevant for you shielding calculation) would clearly speed-up the simulation, e.g. result of the PROFILE card obtained running 1e4 primary with EMFCUT set at 100 keV for e-/e+ and 33 keV for gammas:
# Region Total(ms) % p n e+- g
0 -total- 2.36 100.00 33.0 15.2 10.6 21.3
4 L1BPE 1.05 44.49 37.8 8.3 15.8 6.5
19 TarFront 0.46 19.65 77.1 0.3 15.3 6.9
3 L1Pb 0.45 19.24 0.8 24.4 0.5 63.0
If for your specific case you don’t want to change the EMF thresholds, then I would suggest to repeat the 2-steps approach with the use of importance biasing applied to the different shielding regions (which are already well-suited for this kind of manual biasing: assign an increasing importance to the different layers of the shielding from inside to outside). In this way you can also check if the high Maximum CPU time is actually related to the use of AUTOIMBS.
Finally, as last resort you will need to increase the number of particle simulated.
Oh wow, I have completly missed the PROFILE card so far, but this information is really, really helpful for the optimization process. Thank you very much!
Thanks also for your advice on the next steps and possible approaches. I will try both optimizations with a few primaries and see from the results which is better suited.
