Improve initial state for the weighting potential simulation

[RitaF00]

This PR introduces a new handling of the initial state application to the weighting potential calculation. (#530 )
It fixes the problem of the incorrect p+ weighting potential convergence observed in an undepleted detector when using small values of max_tick_distance (e.g.below 0.3 mm ) to build the grid; the same issue also appeared when passing to the weighting potential calculation a very refined custom grid, such as the grid obtained from the electric potential calculation.

The weighting potential calculated on grids with small max_tick_distance converges too fast; this behaviour is not related to the RCC layer model, since this was also spotted for simulations not including the RCC layer model.

See the example below:
P+ weighting potential for max_tick_distance = 0.5 mm (left), showing correct behaviour, and max_tick_distance = 0.1 mm (right), showing incorrect behaviour.

→ vacuum

→ cryostat

→ RCC layer

To verify if the convergence is correct or not, the minimum of the weighting potential is evaluated in a test volume inside the undepleted region: since this region is in direct contact with the p+ electrode and behaving like a conductive material (undepleted → no charge density), the weighting potential must be equal to 1.
We can see how decreasing the value of max_tick_distance ( i.e. having a denser grid) leads to smaller and smaller weighting potential values.

The new implementation changes the initial state applied for the weighting potential simulation. In the original code, the simulation begins with

apply_initial_state!(sim, potential_type, contact_id,grid)

In this PR, instead of applying the initial state directly to the user grid ( custom grid or build from custom max_tick_distance), the initial state is first applied to the default grid.
The simulation is therefore initialised on a coarse grid and an initial convergence is computed.
Allowing the simulation to perform this step helps the SOR algorithm to reach the correct convergence; when the initial grid is too dense, the algorithm gets stuck around an incorrect solution, and continues to oscillate around it until the iteration stops when maximum of the residual of the weighting potential is smaller than a certain convergence limit ( default 10^(-7)).

In principle, this convergence problem can be solved by forcing the algorithm to iterate for more steps before checking convergence; this can be achieved by decreasing the convergence limit (requiring a smaller error)

or by increasing the number of iterations between checks ( by default fixed to 1000).
Above n_iteration_between_checks = 30.000, the convergence becomes correct.

Naturally, we are not only interested in the correct convergence, but also in achieving it within computationally feasible time frames for our algorithm.

That’s why it’s needed to change how the initial state is handled.

Starting from a coarse grid, helps the convergence to not get stuck and reach the correct convergence. This has been checked using both small max_tick_distance:

and custom grids:

We can now compare the minimum of the weighting potential with the ‘OLD’ and the ‘NEW’ implementation:

We can see that now the minimum of the weighting potential is no longer decreasing, and oscillates around 1; for all the max_tick_distances, the deviation from 1 is always smaller than 0.5%.

The added code performs the following steps (only in the weighting potential simulation):

• Applies the initial state to the coarse default grid
• Updates until convergence
• Map the weighting potential obtained onto the user grid
• Updates until convergence again
• After this, the refinement on the potential is applied as usual.

Introducing this additional convergence step does not add any significant computational time.

[RitaF00]

All tests have been run locally and they passed.

[fhagemann]

Don't worry about the failing tests now, they seem to be related with a new release of Geant4 binaries that require some adjustments in other binary packages.

I combined your four commits into one --> make sure that you pull those changes to your local branch before push here, or you'll might weird merge conflicts.

[fhagemann]

I pushed one more commit (adding comments for future reference, and reducing the point complexity by shortening how we pass keyword arguments).
As the Github tests are broken, let me perform some local tests to make sure this works as intended.
But this looks like a pretty straight-forward change to me!

Thanks for the work and your very detailed description that came with it! 😄

[fhagemann]

@RitaF00 One question I would have:
from your experience, would you recommend ALWAYS doing this, or only when depletion_handling is set to true?

[RitaF00]

@fhagemann From what I’ve seen, the incorrect convergence of the weighting potential only appeared in the case of undepleted detectors on the p+ electrode. Given the detector geometry, this suggests that the issue is mainly related to regions where the potential (or its gradient) changes very rapidly.

In the end, though, introducing this potential initialization does not add any computational overhead; it simply helps “warm up” the SOR algorithm solution.

So in general, I think it does not make much difference whether we apply this method always or only when depletion_handling = true.

Applying it in an agnostic way might actually be the more conservative choice.

What do you think?

[fhagemann]

That’s also what I was thinking, thanks!

Final question: let’s assume the user chooses use the default grid in the calculation: here you would just run the update algorithm twice on the default grid instead of once (essentially doubling the number of iterations on the default grid)?

im just thinking of cases in which the default grid might already have a large number of grid points (e.g. in heavily pixelated detectord), how this might affect the run time

[RitaF00]

mh right.... never worked wiht heavily pixelated detector

In this detectors, is the default grid built as the default case? ( max_tick_distance = axis_length / 4)

In the end we can just apply this 'warm up' only if a custom/user grid is passed.

[fhagemann]

I like this idea of doing the default-grid "warm up" only if a custom grid is passed!

[fhagemann]

The "warm-up" does not seem to care about the max_dick_distance and max_distance_ratio that is passed to the function. Is this intended? If so, we can make a direct comparison between the grid and default_grid and only apply the "warm-up" if they're different.

[fhagemann]

I pushed some changes in how we perform our tests to main that hopefully fix the failing tests due to Geant4 not yet being compatible with the newest Geant4_jll. We'll see if that worked ^^

[fhagemann]

Let's give it a second try..