Conversation
|
I need to see when I can find enough time to look at this in detail. |
|
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 |
|
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 |
src/Simulation/Simulation.jl
Outdated
| electric_potential::Union{ElectricPotential{T}, Missing} | ||
| weighting_potentials::Vector{Any} | ||
| electric_field::Union{ElectricField{T}, Missing} | ||
| spacing_surface_refinement::Union{NTuple{3,T}, Nothing} |
There was a problem hiding this comment.
I'm not a fan of adding this as a field to Simulation.
Could this be something that could be added to the World in sim.world?
There was a problem hiding this comment.
I implemented this
| boundaries: | ||
| left: inf | ||
| right: inf | ||
| spacing_surface_refinement: [1e-4, 1e-4, 1e-4] #m |
There was a problem hiding this comment.
If it's here semantically, I would also add this to the World that we create from it and save it to sim.world.
Not necessarily for this PR, but we could think about also allowing to put things like the refinement limits, max/min tick distance etc. in here, so that they do not have to be passed at runtime to the calculate_potential! methods.
| # Cylindrical: inner loop over z (3rd dimension) | ||
| i1_safe = clamp(i1, 1, size(ϵ_r, 3)) | ||
| in1_safe = clamp(in1, 1, size(ϵ_r, 3)) | ||
| i2_safe = clamp(i2, 1, size(ϵ_r, 2)) | ||
| in2_safe = clamp(in2, 1, size(ϵ_r, 2)) | ||
| i3_safe = clamp(i3, 1, size(ϵ_r, 1)) | ||
| in3_safe = clamp(in3, 1, size(ϵ_r, 1)) | ||
| # ϵ_r is not transformed into an red-black-4D-array. | ||
| # The inner loop (over i1) is along the z-Dimension (Cylindrical Case), | ||
| # which is the 3rd dimension for Cylindrical coordinates: (r, φ, z) | ||
| return @inbounds begin | ||
| ϵ_r[ in3, in2, in1 ], | ||
| ϵ_r[ i3, in2, in1 ], | ||
| ϵ_r[ in3, i2, in1 ], | ||
| ϵ_r[ i3, i2, in1 ], | ||
| ϵ_r[ in3, in2, i1 ], | ||
| ϵ_r[ i3, in2, i1 ], | ||
| ϵ_r[ in3, i2, i1 ], | ||
| ϵ_r[ i3, i2, i1 ] | ||
| ϵ_r[in3_safe, in2_safe, in1_safe], | ||
| ϵ_r[i3_safe, in2_safe, in1_safe], | ||
| ϵ_r[in3_safe, i2_safe, in1_safe], | ||
| ϵ_r[i3_safe, i2_safe, in1_safe], | ||
| ϵ_r[in3_safe, in2_safe, i1_safe], | ||
| ϵ_r[i3_safe, in2_safe, i1_safe], | ||
| ϵ_r[in3_safe, i2_safe, i1_safe], | ||
| ϵ_r[i3_safe, i2_safe, i1_safe] |
There was a problem hiding this comment.
Why is this needed now? Was it needed before already (is this a bug fix)?
Or is this something required now with the things you have added in this PR?
There was a problem hiding this comment.
This bit of the code was giving a bound error when adding more grid points around the surface layer. When refining the grid near the surface, some index would fall just outside the valid array range. Previously this did not show up because the grid configurations we used did not hit these boundary cases, but with the refinement introduced in this PR it can occur and leads to a BoundsError. So the clamping would act as a safeguard to keep the indices within the range. In that sense this is effectively a bug fix for boundary handling, but it only became visible with the refinement changes added in this PR
There was a problem hiding this comment.
Would it make sense to work with the extended grid then, to also keep boundary conditions fulfilled?
SolidStateDetectors.jl/src/Axes/DiscreteAxis.jl
Lines 92 to 123 in f279b07
There was a problem hiding this comment.
This looks like a better approach. I’ll check whether we can use the extended grid here instead
There was a problem hiding this comment.
get_extended_ticks is already used in ϵ_r and all the approaches I tried to fix this issue didn’t work but clamp only affects cases where the index would go out of bounds, so for example the very refined grid of 1e-5m, for all usual cases that previously ran without errors, it changes nothing, so I think it’s a safe fix for very fine grids or edge conditions
There was a problem hiding this comment.
Just by looking at the code, I don't fully understand this?
Can you clarify in the most simple terms why we need this clamping? When refining surface intervals close to the boundaries of the world? If so, why don't we care about the boundary conditions applied to the boundaries of the world?
| # Cartesian: inner loop over x (1st dimension) | ||
| i1_safe = clamp(i1, 1, size(ϵ_r, 1)) | ||
| in1_safe = clamp(in1, 1, size(ϵ_r, 1)) | ||
| i2_safe = clamp(i2, 1, size(ϵ_r, 2)) | ||
| in2_safe = clamp(in2, 1, size(ϵ_r, 2)) | ||
| i3_safe = clamp(i3, 1, size(ϵ_r, 3)) | ||
| in3_safe = clamp(in3, 1, size(ϵ_r, 3)) | ||
| # The inner loop (over i1) is along the x-Dimension (Cartesian Case), | ||
| # which is the 1rd dimension for Cartesian coordinates: (x, y, z) | ||
| return @inbounds begin | ||
| ϵ_r[ in1, in2, in3 ], | ||
| ϵ_r[ in1, in2, i3 ], | ||
| ϵ_r[ in1, i2, in3 ], | ||
| ϵ_r[ in1, i2, i3 ], | ||
| ϵ_r[ i1, in2, in3 ], | ||
| ϵ_r[ i1, in2, i3 ], | ||
| ϵ_r[ i1, i2, in3 ], | ||
| ϵ_r[ i1, i2, i3 ] | ||
| ϵ_r[in1_safe, in2_safe, in3_safe], | ||
| ϵ_r[in1_safe, in2_safe, i3_safe], | ||
| ϵ_r[in1_safe, i2_safe, in3_safe], | ||
| ϵ_r[in1_safe, i2_safe, i3_safe], | ||
| ϵ_r[ i1_safe, in2_safe, in3_safe], | ||
| ϵ_r[ i1_safe, in2_safe, i3_safe], | ||
| ϵ_r[ i1_safe, i2_safe, in3_safe], | ||
| ϵ_r[ i1_safe, i2_safe, i3_safe] |
There was a problem hiding this comment.
I still need to implement this
| for d in 1:3 | ||
| if length(old_grid.axes[d].ticks) != 1 | ||
| if min_spacing[d] < T(1e-5) || min_spacing[d] > T(1e-3) | ||
| @warn "min_spacing[$d] = $(min_spacing[d]) m is out of bounds (1e-5 m < min_spacing < 1e-3 m). Higher values don't require surface refinement" | ||
| if min_spacing[d] < T(1e-5) | ||
| @warn "Using min_spacing[$d] = 1e-5 m" | ||
| else | ||
| @warn "Using min_spacing[$d] = 1e-3 m" | ||
| end | ||
| end | ||
| end | ||
| end |
There was a problem hiding this comment.
How about avoiding throwing a warning for the phi-axis in cylindrical coordinates.
src/Simulation/Simulation.jl
Outdated
| only2d = size(closed_pot.data, 2) == 1 | ||
| int = interpolate_closed_potential(closed_pot, Val(only2d)) |
There was a problem hiding this comment.
Would this perform correctly in Cartesian coordinates?
Would only2d also require some is_cyl check?
There was a problem hiding this comment.
Yes, the 2D interpolation path is only intended for the cylindrical case when the φ dimension has size 1, so maybe it makes more sense to have only2d = is_cyl && size(closed_pot.data, 2) == 1, I added it
src/Simulation/Simulation.jl
Outdated
| pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data, | ||
| not_only_paint_contacts = not_only_paint_contacts, paint_contacts = paint_contacts) |
There was a problem hiding this comment.
You could shorten this to
| pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data, | |
| not_only_paint_contacts = not_only_paint_contacts, paint_contacts = paint_contacts) | |
| pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data; | |
| not_only_paint_contacts, paint_contacts) |
src/Simulation/Simulation.jl
Outdated
| has_surface_model = isdefined(sim.detector.semiconductor.impurity_density_model, :surface_imp_model) | ||
|
|
||
| if has_surface_model | ||
| if sim.spacing_surface_refinement === nothing |
There was a problem hiding this comment.
I usually use isnothing to check equality,
not sure if there is a preferred way:
| if sim.spacing_surface_refinement === nothing | |
| if isnothing(sim.spacing_surface_refinement) |
src/Simulation/Simulation.jl
Outdated
| paint_contacts = paint_contacts, | ||
| sor_consts = is_last_ref ? T(1) : sor_consts ) | ||
|
|
||
| if has_surface_model && iref==3 && sim.spacing_surface_refinement !== nothing |
There was a problem hiding this comment.
| if has_surface_model && iref==3 && sim.spacing_surface_refinement !== nothing | |
| if has_surface_model && iref == 3 && !isnothing(sim.spacing_surface_refinement) |
|
Since there are many changes suggested I commented the ones I've added in the new commits to keep track of them |
|
I think I've now implemented all the changes, let me know if there's still something missing/ or to improve |
|
|
||
| """ | ||
| refine_surface!(sim::Simulation{T}; | ||
| min_spacing::NTuple{3,T} = (T(1e-4), T(1e-4), T(1e-4)), |
There was a problem hiding this comment.
On my end, it still looks like this is not done (?)
| 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)); |
| """ | ||
| struct World{T <: SSDFloat, N, S} <: AbstractWorld{T, N} | ||
| intervals::NTuple{N, SSDInterval{T}} | ||
| spacing_surface_refinement::Union{NTuple{N,T}, Nothing} |
There was a problem hiding this comment.
I think it would be nice to have this type-stable.
Could we replace nothing by (NaN, NaN, NaN), such that it's always an NTuple{N,T}?
| spacing_surface_refinement = if haskey(dict, "spacing_surface_refinement") | ||
| ntuple(i -> _parse_value(T, dict["spacing_surface_refinement"][i], internal_length_unit), 3) | ||
| else | ||
| nothing | ||
| end |
There was a problem hiding this comment.
How do we deal with the fact that the second value for a cylindrical world should not be a "length"? Also: replace nothing by Tuple of NaNs?
| spacing_surface_refinement = if haskey(dict, "spacing_surface_refinement") | |
| ntuple(i -> _parse_value(T, dict["spacing_surface_refinement"][i], internal_length_unit), 3) | |
| else | |
| nothing | |
| end | |
| spacing_surface_refinement = if haskey(dict, "spacing_surface_refinement") | |
| ntuple(i -> _parse_value(T, dict["spacing_surface_refinement"][i], internal_length_unit), 3) | |
| else | |
| ntuple(i -> NaN, 3) | |
| end |
|
|
||
|
|
||
| function CylindricalWorld(T, dict::AbstractDict, input_units::NamedTuple)::World | ||
| function CylindricalWorld(T, dict::AbstractDict, input_units::NamedTuple, spacing_surface_refinement)::World |
There was a problem hiding this comment.
Here (and elsewhere): can you define the type of spacing_surface_refinement?
|
|
||
| old_grid = sim.electric_potential.grid | ||
|
|
||
| # Enforce minimum and maximum spacing for surface refinement (10 µm < min_spacing < 1 mm) |
There was a problem hiding this comment.
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.
| for d in 1:3 | ||
| if length(old_grid.axes[d].ticks) != 1 | ||
| if min_spacing[d] < T(1e-5) || min_spacing[d] > T(1e-3) | ||
| @warn "min_spacing[$d] = $(min_spacing[d]) m is out of bounds (1e-5 m < min_spacing < 1e-3 m). Higher values don't require surface refinement" |
| # Cylindrical: inner loop over z (3rd dimension) | ||
| i1_safe = clamp(i1, 1, size(ϵ_r, 3)) | ||
| in1_safe = clamp(in1, 1, size(ϵ_r, 3)) | ||
| i2_safe = clamp(i2, 1, size(ϵ_r, 2)) | ||
| in2_safe = clamp(in2, 1, size(ϵ_r, 2)) | ||
| i3_safe = clamp(i3, 1, size(ϵ_r, 1)) | ||
| in3_safe = clamp(in3, 1, size(ϵ_r, 1)) | ||
| # ϵ_r is not transformed into an red-black-4D-array. | ||
| # The inner loop (over i1) is along the z-Dimension (Cylindrical Case), | ||
| # which is the 3rd dimension for Cylindrical coordinates: (r, φ, z) | ||
| return @inbounds begin | ||
| ϵ_r[ in3, in2, in1 ], | ||
| ϵ_r[ i3, in2, in1 ], | ||
| ϵ_r[ in3, i2, in1 ], | ||
| ϵ_r[ i3, i2, in1 ], | ||
| ϵ_r[ in3, in2, i1 ], | ||
| ϵ_r[ i3, in2, i1 ], | ||
| ϵ_r[ in3, i2, i1 ], | ||
| ϵ_r[ i3, i2, i1 ] | ||
| ϵ_r[in3_safe, in2_safe, in1_safe], | ||
| ϵ_r[i3_safe, in2_safe, in1_safe], | ||
| ϵ_r[in3_safe, i2_safe, in1_safe], | ||
| ϵ_r[i3_safe, i2_safe, in1_safe], | ||
| ϵ_r[in3_safe, in2_safe, i1_safe], | ||
| ϵ_r[i3_safe, in2_safe, i1_safe], | ||
| ϵ_r[in3_safe, i2_safe, i1_safe], | ||
| ϵ_r[i3_safe, i2_safe, i1_safe] |
There was a problem hiding this comment.
Just by looking at the code, I don't fully understand this?
Can you clarify in the most simple terms why we need this clamping? When refining surface intervals close to the boundaries of the world? If so, why don't we care about the boundary conditions applied to the boundaries of the world?



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:
---> Without RCC layer but using custom grid: same issue was observed:
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: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:
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.