JuliaPhysics/SolidStateDetectors.jl

Grid refinement fails if potential is zero everywhere

fhagemann opened this issue · 0 comments

To test some diffusion and self-repulsion effects (see #209), I tried to define a detector with potential 0V on all contacts.
However, the grid refinement doesn't seem to like this at all.

MWE:

using SolidStateDetectors

sim = Simulation(SSD_examples[:InvertedCoax])
calculate_electric_potential!(sim) # works

sim.detector = SolidStateDetector(sim.detector, contact_id = 2, contact_potential = 0)
calculate_electric_potential!(sim) # fails

throws the following error:

Stacktrace:
  [1] trunc
    @ ./float.jl:716 [inlined]
  [2] floor
    @ ./float.jl:294 [inlined]
  [3] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:631 [inlined]
  [5] getindex
    @ ./broadcast.jl:575 [inlined]
  [6] macro expansion
    @ ./broadcast.jl:984 [inlined]
  [7] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [8] copyto!
    @ ./broadcast.jl:983 [inlined]
  [9] copyto!
    @ ./broadcast.jl:936 [inlined]
 [10] copy
    @ ./broadcast.jl:908 [inlined]
 [11] materialize
    @ ./broadcast.jl:883 [inlined]
 [12] _create_refined_grid(p::ElectricPotential{Float32, 3, Cylindrical, Tuple{SolidStateDetectors.DiscreteAxis{Float32, :r0, :infinite, IntervalSets.ClosedInterval{Float32}}, SolidStateDetectors.DiscreteAxis{Float32, :reflecting, :reflecting, IntervalSets.ClosedInterval{Float32}}, SolidStateDetectors.DiscreteAxis{Float32, :infinite, :infinite, IntervalSets.ClosedInterval{Float32}}}}, max_diffs::Tuple{Float32, Float32, Float32}, minimum_distances::Tuple{Float32, Float32, Float32})
    @ SolidStateDetectors ~/Software/SolidStateDetectors.jl/src/PotentialSimulation/ConvergenceAndRefinement.jl:223
 [13] refine_scalar_potential(p::ElectricPotential{Float32, 3, Cylindrical, Tuple{SolidStateDetectors.DiscreteAxis{Float32, :r0, :infinite, IntervalSets.ClosedInterval{Float32}}, SolidStateDetectors.DiscreteAxis{Float32, :reflecting, :reflecting, IntervalSets.ClosedInterval{Float32}}, SolidStateDetectors.DiscreteAxis{Float32, :infinite, :infinite, IntervalSets.ClosedInterval{Float32}}}}, max_diffs::Tuple{Float32, Float32, Float32}, minimum_distances::Tuple{Float32, Float32, Float32}; only2d::Val{true})
    @ SolidStateDetectors ~/Software/SolidStateDetectors.jl/src/PotentialSimulation/ConvergenceAndRefinement.jl:107
 [14] refine_scalar_potential(p::ElectricPotential{Float32, 3, Cylindrical, Tuple{SolidStateDetectors.DiscreteAxis{Float32, :r0, :infinite, IntervalSets.ClosedInterval{Float32}}, SolidStateDetectors.DiscreteAxis{Float32, :reflecting, :reflecting, IntervalSets.ClosedInterval{Float32}}, SolidStateDetectors.DiscreteAxis{Float32, :infinite, :infinite, IntervalSets.ClosedInterval{Float32}}}}, max_diffs::Tuple{Float32, Float32, Float32}, minimum_distances::Tuple{Float32, Float32, Float32})
    @ SolidStateDetectors ~/Software/SolidStateDetectors.jl/src/PotentialSimulation/ConvergenceAndRefinement.jl:106
 [15] refine!(sim::Simulation{Float32, Cylindrical}, ::Type{ElectricPotential}, max_diffs::Tuple{Float32, Float32, Float32}, minimum_distances::Tuple{Float32, Float32, Float32}; not_only_paint_contacts::Bool, paint_contacts::Bool, update_other_fields::Bool)
    @ SolidStateDetectors ~/Software/SolidStateDetectors.jl/src/Simulation/Simulation.jl:703
 [16] refine!
    @ ~/Software/SolidStateDetectors.jl/src/Simulation/Simulation.jl:703 [inlined]
 [17] _calculate_potential!(sim::Simulation{Float32, Cylindrical}, potential_type::UnionAll, contact_id::Missing; grid::Missing, convergence_limit::Float64, refinement_limits::Vector{Float64}, min_tick_distance::Missing, max_tick_distance::Missing, max_distance_ratio::Int64, depletion_handling::Bool, use_nthreads::Int64, sor_consts::Missing, max_n_iterations::Int64, n_iterations_between_checks::Int64, not_only_paint_contacts::Bool, paint_contacts::Bool, verbose::Bool)
    @ SolidStateDetectors ~/Software/SolidStateDetectors.jl/src/Simulation/Simulation.jl:884
 [18] _calculate_potential!(sim::Simulation{Float32, Cylindrical}, potential_type::UnionAll, contact_id::Missing)
    @ SolidStateDetectors ~/Software/SolidStateDetectors.jl/src/Simulation/Simulation.jl:768
 [19] #calculate_electric_potential!#165
    @ ~/Software/SolidStateDetectors.jl/src/Simulation/Simulation.jl:1054 [inlined]
 [20] calculate_electric_potential!(::Simulation{Float32, Cylindrical})
    @ SolidStateDetectors ~/Software/SolidStateDetectors.jl/src/Simulation/Simulation.jl:1054
 [21] top-level scope
    @ In[8]:1
 [22] eval
    @ ./boot.jl:360 [inlined]
 [23] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1090

I assume that this comes from the first lines of _create_refined_grid:

function _create_refined_grid(p::ScalarPotential{T,3}, max_diffs::NTuple{3, T}, minimum_distances::NTuple{3, T}) where {T}
    n_1 = floor.(Int, [maximum(abs.(p.data[i+1,:,:] .- p.data[i,:,:])) for i in 1:size(p.data, 1)-1] ./ max_diffs[1]) 
    n_2 = floor.(Int, [maximum(abs.(p.data[:,i+1,:] .- p.data[:,i,:])) for i in 1:size(p.data, 2)-1] ./ max_diffs[2]) 
    n_3 = floor.(Int, [maximum(abs.(p.data[:,:,i+1] .- p.data[:,:,i])) for i in 1:size(p.data, 3)-1] ./ max_diffs[3]) 
    ns = (n_1, n_2, n_3)
    #...
end

where max_diffs = (0, 0, 0) if the sim.electric_potential is just 0 everywhere, resulting in NaN for n_1, n_2 and n_3.
In addition, we should also avoid max_diffs default being 0 in refine!:

function refine!(sim::Simulation{T}, ::Type{WeightingPotential}, contact_id::Int,

                    max_diffs::Tuple{<:Real,<:Real,<:Real} = (T(0), T(0), T(0)),  # should be avoided !!!

                    minimum_distances::Tuple{<:Real,<:Real,<:Real} = (T(0), T(0), T(0))) where {T <: SSDFloat}
    sim.weighting_potentials[contact_id] = refine_scalar_potential(sim.weighting_potentials[contact_id], max_diffs, minimum_distances)
    nothing
end

We should make sure that max_diffs is never (0,0,0) or that the grid is just not refined when this is the case.