I’ve recently pulled the most recent Garfield++ version from the gitlab repository and I noticed that my DriftLineRKF simulations using COMSOL field maps are now taking noticeably longer. An individual event that would have taken a few seconds to simulate now takes more than 10 minutes (and sometimes much more). I’ve also been getting a lot of DriftLineRKF::IntegrateAlpha and DriftLineRKF::Avalanche warnings and a lot more mesh convergence warnings, none of which were happening before the update.
So far, I have tried:
Updating my .gas files (from v12 to v13);
Disabling the signal calculation;
Updating ROOT to the latest stable release (6.32.12);
Removing the #pragma omp parallel for directive from the DriftLineRKF loop;
Simplifying the COMSOL model and mesh;
The simpler COMSOL model/mesh did help with the DriftLineRKF warnings but that’s about it. It still takes a large amount of time to run a single event. I also had to remove from the model the electrodes used for the signal calculation, which isn’t ideal.
I rebuilt Garfield++ using the old hash and everything went back to normal. The hash of the version is c9440b06fc (merge ‘gpu-signal-calc’, Aug. 05, 2024). The events in this version should finish pretty quickly with OpenMP, and with very few warnings too (or none at all depending on the maximum step size).
I’ve attached a reproducible version of the problem. The COMSOL field map/mesh are from the simpler version, which is a single zero-centered, 20[um], 1400[V] wire inside a box of air with grounded surfaces.
Some good news: I looked at the commit history of the repo again and found the merge of a branch named “COMSOLFieldMap” into master. The commit before that merge (“Default set to importing all materials”, SHA 3891476306) could be the one causing the behavior I described. Changing the head to b4757d7ba1 (“m_GridMesh → m_grid”) retains the behavior of the older version I was using. I’ll check again in the morning.
Hello, here is some additional information I gathered:
This is a side by side comparison of the last std output of the original simulation I’m working on (includes signal calculations, weighting fields and ROOT I/O classes). On the left terminal, the HEAD of the Garfield++ repo is at 3891476306 (“m_GridMesh → m_grid”), while the right terminal used the most recent version. I ran the exact same piece of code on both terminals:
The “Average time per event” output includes the time it takes to read the COMSOL files (2.8GB+, fields and mesh), so the actual drift simulation on the older version is, on average, faster than 27s. Moreover, the results returned by the functions DriftLineRKF::GetAvalancheSize() and DriftLineRKF::GetLoss() are returning a large amount of NaN in the most recent version. Take for example the number of electrons returned by DriftLineRKF::GetAvalancheSize() in one of the 10 events I simulated, it has 223 empty entries:
y = Interpolate1D(e, table[0][0], m_eFields, intp, extr, logval);
to
y = Interpolate1D(e, table[0][0], m_eFields, intp, extr, false);
and check if you get back the “old” results?
The commit changes the way the Townsend coefficient and attachment coefficient are extrapolated towards high (and low) electric fields. If this has a significant impact on your results, it suggests that the range of electric fields in the gas file you are is not large enough to cover the range of electric fields in your detector.
Thank you so much! I hadn’t even considered the electric field range in the .gas file. I went to check the COMSOL model for the largest magnitude of the electric field as soon as I finished reading your post. Since the wire is so thin, some of the regions in the detector have fields of 500+ kV/cm, while the gas file I was using only covered the [0.1, 100] kV/cm range (embarrassing! haha). Using a file that covered the [0.1, 1000] kV/cm range did the trick and, as one would expect, improved the results all around.
If this is not the expected behavior, then there seems to be some sort of problem with the exponential extrapolation of the attachment coefficient, which would explain my previous problem with DriftLineRKF::GetLoss() and DriftLineRKF::GetAvalancheSize(). I’ve attached at the end of the post the electron attachment coefficient table exported by ViewMedium::EnableExport() as well as the gas file.
If you take a look at the end of the electronEta.txt file you’ll see a few NaN for the extrapolated fields:
The problem seems to be consistent with the extrapolation towards larger fields, since calculating a new gas table that covered a larger range solved my previous problem.
If I understood the code correctly, the exponential extrapolation is calculated in the lines 1503-1508 of Medium.cc (since m_extrAtt.second == 1, line 748 in Medium.hh):
else {
// Log values in gas table for alpha, eta.
const double extr2 = (std::exp(ytab[nt - 1]) - std::exp(ytab[nt - 2])) /
(xtab[nt - 1] - xtab[nt - 2]);
result = log(std::exp(ytab[nt - 1]) + extr2 * (x - xtab[nt - 1]));
}
This means that we’ll have extr2 < 0 for every ytab[nt - 1] < ytab[nt - 2]. That means that for very large values of x the program will have to compute the natural logarithm of a negative value since abs(extr2 * (x - xtab[nt - 1])) > exp(ytab[nt - 1]).
So the problem can be summed up to ytab[nt - 1] < ytab[nt - 2], logval = true and x >>> xtab[n -1]
So if we’re looking at the eta in the gas file “ar_97_co2_3.gas” from the Examples/GasFile folder:
The graph shows that for x = 190.7E3 the value of y is undefined, which matches the output of the modified version of the printTable.C example I just ran:
Thanks @gabrielribcesario@hschindl for spotting this issue. It indeed makes sense to avoid to calculate log(0). The log calculation is likely expensive, so I would perform a check before.
I guess these calculations should somehow be optimized, especially for very low values where it does not make sense to spend a lot of computing time to calculate a number that is close to zero. We should maybe sit down and think whether to keep the alpha and eta in logscale in the gastables.
Apologies for the discomfort
greets
Piet
