Surface Grid Refinement

[claudiaalvgar]

This PR introduces surface-only grid refinement for electric potential simulations and fixes boundary-related indexing issues that could lead to segmentation faults at high refinement resolutions e.g for 1e-5m grid spacing. In addition, it also fixes an existent issue when calculating the electric potential with custom grids for small grid spacing. For custom small grid spacing (e.g. 2e-5m), the detector appeared to be undepleted while for bigger grid spacing (1e-4m) the behaviour was correct. This issue was not related to the new RCC layer model since this was also spotted for simulations not including the RCC layer model. See the issue below:
---> With RCC layer custom initial grid:

• 0.02mm grid spacing: wrong behaviour (I attach the code in case this issue has to be reproduced):

T = Float64 
sim = Simulation{T}(SSD_examples[:IVCIlayer])
calculate_electric_potential!(sim ,max_n_iterations = 10 ,grid= Grid(sim) ,verbose = true ,depletion_handling = true)
g = sim.electric_potential.grid
ax1, ax2, ax3 = g.axes
tick_dis=0.02*mm 
user_additional_ticks_ax1=sort(vcat(ax1.interval.left:tick_dis:ax1.interval.right)) 
user_additional_ticks_ax3=sort(vcat(ax3.interval.left:tick_dis:ax3.interval.right))[2:end] # only support even number of ticks in z-direction 
user_ax1 = typeof(ax1)(ax1.interval, user_additional_ticks_ax1) 
user_ax3 = typeof(ax3)(ax3.interval, user_additional_ticks_ax3) 
user_g = typeof(g)((user_ax1, ax2, user_ax3)) 
calculate_electric_potential!(sim, refinement_limits = [0.2, 0.1, 0.05, 0.02], use_nthreads=8, grid = user_g, depletion_handling = true)

* 0.1mm grid spacing custom initial grid: correct behaviour:

---> Without RCC layer but using custom grid: same issue was observed:

• 0.02mm custom grid spacing: wrong behaviour

* 0.1mm custom grid spacing: correct behaviour The new implementation checks whether a given slice of point types contains any surface (inactive layer) points and builds a new refined grid across all three axes. Surface refinement is automatically triggered during the electric potential calculation when:

• A surface impurity model is present

• At least 3 refinement steps are used in the electric potential calculation

• Final refinement limit ≤ 0.05

A YAML entry was added to configure spacing_surface_refinement, that controls the biggest grid spacing allowed in the surface for the 3 axis. The new method identifies the surface regions and it just refines in this area as can be seen, for example for a 0.1mm spacing:
Having this spacing_surface_refinement: [1e-4, 1e-4, 1e-4] in the yaml file and running the default: calculate_electric_potential!(sim, use_nthreads=8, grid=g, verbose = true, depletion_handling = true) will now give the correct behaviour for the surface layer, or even: calculate_electric_potential!(sim, refinement_limits = [0.2, 0.1, 0.05, 0.02], use_nthreads=8, grid = user_g, depletion_handling = true) if we want to further refine the bulk:

The weighting potential is also well behaved even when passing the custom grid from the surface-refined electric potential:

This provides a faster calculation O(3min) than for the custom grid case O(10min) since the initial grid was refined from the very beginning, see old case behaviour for custom grid:

The surface refinement is applied after three initial bulk refinements (the default refinement setup for the electric potential calculation) because only at that stage does the surface behaviour become well-resolved. Therefore, this model requires the default bulk refinement setup in order to correctly apply the surface refinement: Default: refinement_limits = [0.2, 0.1, 0.05]) + surface refinement. Applying the surface refinement after the first or second refinement step does not fully capture the surface layer, leading to incomplete refinement, as you can see :

• Surface Refinement after 2nd refinement uses [0.2, 0.1, -surface ref- 0.05]) and min_spacing = (1e-4, 1e-4, 1e-4)m in surface layer

• Initial surface refinement (this uses -surface ref- + [0.2, 0.1, 0.05]) and min_spacing = (1e-4, 1e-4, 1e-4)m in surface layer

This PR implements the default behaviour (3 refinements) plus the surface refinement, ensuring accurate resolution near surfaces. Additionally, the model now supports refining surface intervals down to 1e-5m, preserving correct physical behaviour and avoiding the incorrect results that can occur with custom grids resulting in undepleted detectors.

• New method 0.02mm spacing: Correct behaviour

NOTE: this model has not been tested on cartesian grid defined detectors but I wanted to ask for inputs before. I wasn't sure neither if this should have been a Draft PR, so feel free to move it.

[fhagemann]

I need to see when I can find enough time to look at this in detail.
Are there any tests you could add to test the new code for it's correct behavior (When it should be throwing warnings, edge cases, etc.) that you could add here?

[claudiaalvgar]

Sure no problem. This code is yet a bit preliminary since it doesn't handle the 3D case i.e. it won't refine the phi axis (this I found a bit tricky) and I also didn't try the model on cartesian detectors but I think this should work

[claudiaalvgar]

Update: when doing some tests with other detectors, I saw that refining only inside the surface was sometimes not enough since it didn't capture the edges very well. See this example for a BEGe detector where the surface goes beyond the detector limit:

This was solved in my last commit by refining also in a margin around the surface layer (the variables extra_before and extra_after were added, that I set to be 5 intervals around the surface), with this, we get a much better definition of the surface:

[claudiaalvgar]

I've now tested this model on cartesian defined detectors and also on partial-phi defined detectors. The surface refinement is disabled in phi angle. I also added tests for edge cases and I think if all tests pass we can move this to normal PR review

[fhagemann]

Same thing here.

[claudiaalvgar]

I still need to implement this

[claudiaalvgar]

Since there are many changes suggested I commented the ones I've added in the new commits to keep track of them

[claudiaalvgar]

I think I've now implemented all the changes, let me know if there's still something missing/ or to improve

[fhagemann]

On my end, it still looks like this is not done (?)

Suggested change

	min_spacing::NTuple{3,T} = (T(1e-4), T(1e-4), T(1e-4)),
	min_spacing::NTuple{3,T} = (T(1e-4), T(1e-4), T(1e-4));

[fhagemann]

I think it would be nice to make this depend on the dimensions of the world.
I wouldn't be surprised if smaller/larger worlds would be better off with smaller/larger values.

[fhagemann]

See previous comment

[claudiaalvgar]

The _safe indices in get_ϵ_of_oktant functions I needed to introduce to prevent out-of-bounds access. For very small surface refinements, the code probes points slightly beyond the limits, so some indices like in1, in2, in3 can fall outside the valid range, I don’t know why the code behaves like this but I tried other solutions and none of them worked.
Clamping just forces those indices back inside the grid so we can safely read ϵᵣ. The issue only appears for small refinements because for the rest of refinements the boundaries are within limits, you can reproduce this issue if you want by running for spacing_surface_refinement: [1e-5, 1e-5, 1e-5] in ivc_inactivelayer.yaml with the default get_ϵ_of_oktant functions (it’ll give a SegFault after approximately an hour since calculating the electric potential is slow for very small grids). I’m open for other solutions I just couldn’t find any other better given the logic of the code

[fhagemann]

The _safe indices in get_ϵ_of_oktant functions I needed to introduce to prevent out-of-bounds access. For very small surface refinements, the code probes points slightly beyond the limits, so some indices like in1, in2, in3 can fall outside the valid range, I don’t know why the code behaves like this but I tried other solutions and none of them worked. Clamping just forces those indices back inside the grid so we can safely read ϵᵣ. The issue only appears for small refinements because for the rest of refinements the boundaries are within limits, you can reproduce this issue if you want by running for spacing_surface_refinement: [1e-5, 1e-5, 1e-5] in ivc_inactivelayer.yaml with the default get_ϵ_of_oktant functions (it’ll give a SegFault after approximately an hour since calculating the electric potential is slow for very small grids). I’m open for other solutions I just couldn’t find any other better given the logic of the code

And this is because, instead of expecting maximum 1 tick going "outside" of the grid, you might have more than 1 now?

[fhagemann]

I just want to make sure that we expect accessing ϵ_r out of bounds of not..
Sure, this might fix the fact that we will not get an error, but maybe we shouldn't exceed the bounds of ϵ_r at all in the first place 🤷‍♂️

[fhagemann]

Right now, there is no possibility to pass extra_before and extra_After to _refine_axis_surface when calling refine_surface!.
Is this intended?

[fhagemann]

Do you have an example (some config file + some simulation settings when running calculate_electric_potential! that would cause the calculation to fail when not including these "safe" clamps for ϵ_r.
I might want to have a short look at that myself to see what's going on :)

[fhagemann]

We tracked down the out-of-bounds error when accessing ϵ to the grid not having an even number of grid ticks.
I recall that either all, or just one axis needs to have an even number of ticks for the red-black division of the grid to work. Having an even number of ticks resolves this issue.

Original code:

Forcing the z-axis to have an even number of ticks:

I can look into implementing a fix for this later today or this week.

[claudiaalvgar]

The new changes look good to me and the code behaves as intended so I think this is good to go on my side

[fhagemann]

Okay, I will deal with min/max distance and also look into performance, but will merge this afterwards!

[fhagemann]

surface_intervals are not used anymore after calling this functions, so we should be fine modifying this array instead of allocating more memory.

[fhagemann]

Instead of creating a merged array of unit ranges, we can merge those two for-loops into one and just perform interval refinement for the i that have refine_flags set.

[fhagemann]

No need to allocate memory for sub_widths if the values can just be calculated "on-the-fly" in the for-loop

[fhagemann]

This avoids allocations by preventing creating a Vector for max_spacing, but letting it remain it a Tuple throughout the function.

[claudiaalvgar]

Final comment: I tested the code with new implementations on CPU and GPU and the grid results are as expected, I leave them here for future reference:

• max_spacing = 1e-4 m IVCIlayer on GPU: Elapsed time: 47.767041562 seconds GRID: Grid{Float64, 3, Cylindrical}(0.0 .. 0.04 - length = 348, 0.0 .. 0.0 - length = 1, -0.01 .. 0.09 - length = 902)

* max_spacing = 2.5e-5 m IVCIlayer on GPU: Elapsed time: 80.683083286 seconds GRID: Grid{Float64, 3, Cylindrical}(0.0 .. 0.04 - length = 1292, 0.0 .. 0.0 - length = 1, -0.01 .. 0.09 - length = 3438) * max_spacing = 1e-5 m IVCIlayer on GPU: Elapsed time: 275.473297317 seconds GRID: Grid{Float64, 3, Cylindrical}(0.0 .. 0.04 - length = 3152, 0.0 .. 0.0 - length = 1, -0.01 .. 0.09 - length = 8520)

• max_spacing = 1e-4 m IVCIlayer on CPU Threads= 45: Elapsed time: 28.721281699 seconds GRID: Grid{Float64, 3, Cylindrical}(0.0 .. 0.04 - length = 348, 0.0 .. 0.0 - length = 1, -0.01 .. 0.09 - length = 902)

• max_spacing = 2.5e-5 m IVCIlayer on CPU Threads= 45: Elapsed time: 190.124053047 seconds GRID: Grid{Float64, 3, Cylindrical}(0.0 .. 0.04 - length = 1292, 0.0 .. 0.0 - length = 1, -0.01 .. 0.09 - length = 3438)

• max_spacing = 1e-5 m IVCIlayer on CPU Threads= 45: Elapsed time: 1213.737185377 seconds GRID: Grid{Float64, 3, Cylindrical}(0.0 .. 0.04 - length = 3152, 0.0 .. 0.0 - length = 1, -0.01 .. 0.09 - length = 8520)