search

JuliaPhysics / SolidStateDetectors.jl

Solid state detector field and charge drift simulation in Julia
Other
82 stars 32 forks source link

`:samplesurface` does not work when important points are not in the final geometry #256

Closed fhagemann closed 2 years ago

fhagemann commented 2 years ago

The function get_scale that scales the distance between points for :samplesurface plotting does not work for geometries, for which the important points are not in the final geometry (excluded in CSGDifference or CSGIntersection).

Minimum working example

Detector configuration file

units:
  length: mm
  angle: deg
grid: 
  coordinates: cartesian
  axes:
    x: 
      from: -3
      to: 3
    y:
      from: -3
      to: 3
    z: 
      from: -3
      to: 3
medium: vacuum
detectors:
- semiconductor:
    material: HPGe
    geometry:
      intersection:
        - box: 
            widths: [2,2,2]
        - sphere:
            r: 1.41421356237
  contacts:
    - id: 1
      potential: 100
      geometry:
        tube:
          r: 1
          h: 0
          origin: [0,0,1]
    - id: 2
      potential: 0
      geometry:
        tube:
          r: 1
          h: 0
          origin: [0,0,-1]

Plotting the detector works for :csg, :mesh3d and :wireframe:

using SolidStateDetectors
using Plots; gr(fmt = :png)
sim = Simulation("MWEsamplesurface.yaml");
plot(sim.detector.semiconductor, size = (500,500), st = :mesh3d, color = :lightgray, lw = 0.5)
plot!(sim.detector)

csg

but fails for :samplesurface:

plot(sim.detector, show_semiconductor = true, st = :samplesurface, size = (500,500), n_samples = 50)
ArgumentError: reducing over an empty collection is not allowed

Stacktrace:
  [1] _empty_reduce_error()
    @ Base ./reduce.jl:301
  [2] reduce_empty(op::Function, #unused#::Type{Float32})
    @ Base ./reduce.jl:311
  [3] mapreduce_empty(#unused#::typeof(identity), op::Function, T::Type)
    @ Base ./reduce.jl:345
  [4] reduce_empty(op::Base.MappingRF{typeof(identity), typeof(max)}, #unused#::Type{Float32})
    @ Base ./reduce.jl:331
  [5] reduce_empty_iter
    @ ./reduce.jl:357 [inlined]
  [6] mapreduce_empty_iter(f::Function, op::Function, itr::Vector{Float32}, ItrEltype::Base.HasEltype)
    @ Base ./reduce.jl:353
  [7] _mapreduce(f::typeof(identity), op::typeof(max), #unused#::IndexLinear, A::Vector{Float32})
    @ Base ./reduce.jl:402
  [8] _mapreduce_dim
    @ ./reducedim.jl:330 [inlined]
  [9] #mapreduce#725
    @ ./reducedim.jl:322 [inlined]
 [10] mapreduce
    @ ./reducedim.jl:322 [inlined]
 [11] #_maximum#743
    @ ./reducedim.jl:894 [inlined]
 [12] _maximum
    @ ./reducedim.jl:894 [inlined]
 [13] #_maximum#742
    @ ./reducedim.jl:893 [inlined]
 [14] _maximum
    @ ./reducedim.jl:893 [inlined]
 [15] #maximum#740
    @ ./reducedim.jl:889 [inlined]
 [16] maximum(a::Vector{Float32})
    @ Base ./reducedim.jl:889
 [17] get_scale(csg::SolidStateDetectors.ConstructiveSolidGeometry.CSGIntersection{Float32, SolidStateDetectors.ConstructiveSolidGeometry.Box{Float32, SolidStateDetectors.ConstructiveSolidGeometry.ClosedPrimitive}, SolidStateDetectors.ConstructiveSolidGeometry.Ellipsoid{Float32, SolidStateDetectors.ConstructiveSolidGeometry.ClosedPrimitive, Float32, Nothing, Nothing}})
    @ SolidStateDetectors.ConstructiveSolidGeometry ~/Software/SolidStateDetectors.jl/src/ConstructiveSolidGeometry/CSG.jl:223
 [18] (::SolidStateDetectors.var"#393#394")(o::SolidStateDetectors.Semiconductor{Float32, SolidStateDetectors.ConstructiveSolidGeometry.CSGIntersection{Float32, SolidStateDetectors.ConstructiveSolidGeometry.Box{Float32, SolidStateDetectors.ConstructiveSolidGeometry.ClosedPrimitive}, SolidStateDetectors.ConstructiveSolidGeometry.Ellipsoid{Float32, SolidStateDetectors.ConstructiveSolidGeometry.ClosedPrimitive, Float32, Nothing, Nothing}}, NamedTuple{(:E_ionisation, :f_fano, :ϵ_r, :ρ, :name, :ml, :mt, :diffusion_fieldvector_electrons, :diffusion_fieldvector_holes), Tuple{Unitful.Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, Float64, Float64, Unitful.Quantity{Float64, 𝐌 𝐋^-3, Unitful.FreeUnits{(g, cm^-3), 𝐌 𝐋^-3, nothing}}, String, Float64, Float64, Int64, Int64}}, ElectricFieldChargeDriftModel{Float32}, SolidStateDetectors.ConstantImpurityDensity{Float32}})
    @ SolidStateDetectors ./none:0
 [19] iterate
    @ ./generator.jl:47 [inlined]
 [20] collect_to!(dest::Vector{Float32}, itr::Base.Generator{Vector{Any}, SolidStateDetectors.var"#393#394"}, offs::Int64, st::Int64)
    @ Base ./array.jl:760
 [21] collect_to_with_first!(dest::Vector{Float32}, v1::Float32, itr::Base.Generator{Vector{Any}, SolidStateDetectors.var"#393#394"}, st::Int64)
    @ Base ./array.jl:738
 [22] collect(itr::Base.Generator{Vector{Any}, SolidStateDetectors.var"#393#394"})
    @ Base ./array.jl:719
 [23] macro expansion
    @ ~/Software/SolidStateDetectors.jl/src/PlotRecipes/SolidStateDetector.jl:64 [inlined]
 [24] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, det::SolidStateDetector)
    @ SolidStateDetectors ~/.julia/packages/RecipesBase/qpxEX/src/RecipesBase.jl:289
 [25] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/7ijBv/src/user_recipe.jl:36
 [26] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/7ijBv/src/RecipesPipeline.jl:70
 [27] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
    @ Plots ~/.julia/packages/Plots/YAlrZ/src/plot.jl:208
 [28] #plot#143
    @ ~/.julia/packages/Plots/YAlrZ/src/plot.jl:91 [inlined]
 [29] top-level scope
    @ In[13]:1
 [30] eval
    @ ./boot.jl:373 [inlined]
 [31] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1196
fhagemann commented 2 years ago

It looks like this is because the filter! method in get_scale results in points being empty. https://github.com/JuliaPhysics/SolidStateDetectors.jl/blob/f1bc27859b38a1e70d61a1f1826964f137dded8d/src/ConstructiveSolidGeometry/CSG.jl#L221

A quick-and-dirty fix would be to add

if isempty(points) points = [point for p in primitives(csg) for s in surfaces(p) for point in vertices(s, 4)] end

in the line after, which results in samplesurface

But there might be more elegant solutions to solve this.

fhagemann commented 2 years ago

There also seems to be a problem with get_world_limits_from_objects if theaxes of the grid are not specified in the configuration file but determined from the objects.

Minimum working example

Detector configuration file:

units:
  length: mm
  angle: deg
grid: cartesian
medium: vacuum
detectors:
- semiconductor:
    material: HPGe
    geometry:
      intersection:
        - box: 
            widths: [2,2,2]
        - sphere:
            r: 1.41421356237
  contacts:
    - id: 1
      potential: 100
      geometry:
        tube:
          r: 1
          h: 0
          origin: [0,0,1]
    - id: 2
      potential: 0
      geometry:
        tube:
          r: 1
          h: 0
          origin: [0,0,-1]
using SolidStateDetectors
sim = Simulation("MWE.yaml")
DimensionMismatch("arrays could not be broadcast to a common size")

Stacktrace:
  [1] _bcs1(a::StaticArrays.SOneTo{4}, b::StaticArrays.SOneTo{6})
    @ StaticArrays ~/.julia/packages/StaticArrays/DkWdL/src/broadcast.jl:44
  [2] _bcs(shape::Tuple{StaticArrays.SOneTo{4}}, newshape::Tuple{StaticArrays.SOneTo{6}})
    @ Base.Broadcast ./broadcast.jl:510
  [3] broadcast_shape(::Tuple{StaticArrays.SOneTo{4}}, ::Tuple{StaticArrays.SOneTo{6}})
    @ Base.Broadcast ./broadcast.jl:504
  [4] combine_axes
    @ ./broadcast.jl:499 [inlined]
  [5] combine_axes (repeats 5 times)
    @ ./broadcast.jl:498 [inlined]
  [6] instantiate
    @ ./broadcast.jl:281 [inlined]
  [7] materialize(bc::Base.Broadcast.Broadcasted{StaticArrays.StaticArrayStyle{1}, Nothing, typeof(hcat), Tuple{StaticArrays.SVector{4, CartesianPoint{Float32}}, StaticArrays.SVector{4, CartesianPoint{Float32}}, StaticArrays.SVector{4, CartesianPoint{Float32}}, StaticArrays.SVector{4, CartesianPoint{Float32}}, StaticArrays.SVector{4, CartesianPoint{Float32}}, StaticArrays.SVector{4, CartesianPoint{Float32}}, StaticArrays.SVector{6, CartesianPoint{Float32}}}})
    @ Base.Broadcast ./broadcast.jl:904
  [8] get_world_limits_from_objects(#unused#::Type{Cartesian}, det::SolidStateDetector{Float32, SolidStateDetectors.Semiconductor{Float32, SolidStateDetectors.ConstructiveSolidGeometry.CSGIntersection{Float32, SolidStateDetectors.ConstructiveSolidGeometry.Box{Float32, SolidStateDetectors.ConstructiveSolidGeometry.ClosedPrimitive}, SolidStateDetectors.ConstructiveSolidGeometry.Ellipsoid{Float32, SolidStateDetectors.ConstructiveSolidGeometry.ClosedPrimitive, Float32, Nothing, Nothing}}, NamedTuple{(:E_ionisation, :f_fano, :ϵ_r, :ρ, :name, :ml, :mt, :diffusion_fieldvector_electrons, :diffusion_fieldvector_holes), Tuple{Unitful.Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, Float64, Float64, Unitful.Quantity{Float64, 𝐌 𝐋^-3, Unitful.FreeUnits{(g, cm^-3), 𝐌 𝐋^-3, nothing}}, String, Float64, Float64, Int64, Int64}}, ElectricFieldChargeDriftModel{Float32}, SolidStateDetectors.ConstantImpurityDensity{Float32}}, Vector{SolidStateDetectors.Contact{Float32, SolidStateDetectors.ConstructiveSolidGeometry.Cone{Float32, SolidStateDetectors.ConstructiveSolidGeometry.ClosedPrimitive, Float32, Nothing}, NamedTuple{(:E_ionisation, :f_fano, :ϵ_r, :ρ, :name, :ml, :mt, :diffusion_fieldvector_electrons, :diffusion_fieldvector_holes), Tuple{Unitful.Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, Float64, Float64, Unitful.Quantity{Float64, 𝐌 𝐋^-3, Unitful.FreeUnits{(g, cm^-3), 𝐌 𝐋^-3, nothing}}, String, Float64, Float64, Int64, Int64}}}}, Missing, Missing})
    @ SolidStateDetectors ~/Software/SolidStateDetectors.jl/src/SolidStateDetector/SolidStateDetector.jl:79
  [9] (Simulation{Float32})(dict::Dict{Any, Any})
    @ SolidStateDetectors ~/Software/SolidStateDetectors.jl/src/Simulation/Simulation.jl:161
 [10] Simulation
    @ ~/Software/SolidStateDetectors.jl/src/Simulation/Simulation.jl:171 [inlined]
 [11] Simulation(config_file::String)
    @ SolidStateDetectors ~/Software/SolidStateDetectors.jl/src/Simulation/Simulation.jl:174
 [12] top-level scope
    @ In[2]:1
 [13] eval
    @ ./boot.jl:373 [inlined]
 [14] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1196