Hello everyone,

I was wondering whether it's possible to add the plotting of the full geometry rather than just the contacts. We want to implement a somewhat more involved geometry and debugging correct positions is quite tricky without visualization or boundary checking methods. I assume defining everything as a contact for debugging would be considered bad practice, but something like that would be of great help.

Thanks!

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

PackageCompiler now crashes when trying to compile a Julia system image that contains SolidStateDetectors and a simple precompile execution file like

using SolidStateDetectors

detector_config_filename = SSD_examples[:InvertedCoax]

# PackageCompiler crashes here:
simulation = Simulation{Float32}(detector_config_filename)

This used to work fine - not sure if this is due to changes in SolidStateDetectors, a change in one of it's depencencies, or Julia v1.4 -> v1.5.

When improving the application of boundary conditions (#206) I recognized that
in(pt::AbstractCoordinatePoint, g::AbstractGeometry) is not allocation free anymore
if the geometry consists out of too many (>4) primitives.

E.g., in case of the geometry defining the mantle contact of the IVC example detector:

using SolidStateDetectors, BenchmarkTools
sim = Simulation(SSD_examples[:InvertedCoax]);
g = sim.detector.contacts[2].geometry; # consists out of 6 primitives
pt = zero(CartesianPoint{Float32});

@btime in($pt, $g); # -> Not allocation free
  148.439 ns (6 allocations: 608 bytes)

For geometries with 4 or less primitives no allocations are caused:

gAllocFree = g.a.a; # consists out of 4 primitives
@btime in($pt, $gAllocFree); 
  58.282 ns (0 allocations: 0 bytes)

As one can see this also worsens the performance.
This causes also some major amount of allocations in the application of boundary conditions (#85).
It will also be an issue in the charge drift where this function is also used quite a lot of times.

I am not sure whether 4 -> 5 is always the edge. It might also depend on the types of the primitives.

In addition, the number of allocations increases for increasing number of primitives of the geometry:

@btime in($pt, $(g.a)); # -> Only 3 allocations
105.006 ns (3 allocations: 272 bytes)

The excluded individual primitives (g -> gAllocFree) are allocations free:

@btime in($pt, $(g.a.b));
  11.071 ns (0 allocations: 0 bytes)

and

@btime in($pt, $(g.b));
  11.371 ns (0 allocations: 0 bytes)

Before the electric potential is calculated via SOR, the initial boundary conditions are applied onto the grid
and some constant weights (needed in the SOR) are precalculated. This also happens after each grid refinement.

If there are many objects in the world, e.g. highly pixelated detectors, and the grid becomes finer and finer,
these precalculations actually can take quite some time.

There are some ways to improve this:

• Get rid of unnecessary point conversions
  The time critical part is the lookup to which object each individual point belongs: the in-method.
  Right now, we only parse the grid point in form of the coordinate system of the world, e.g. cylindrical coordinates, to this method. Depending on the geometry of the object, the parsed point might be converted into the other coordinate system, e.g. cartesian. And this several times. We should do this just once and parse the point in both forms to the in methods.

• Parallize (#206)
  Right now, no multithreading is used here.
  But the loop running over the grid points can be parallized.

• Remove allocations in in method for larger CSG's
  See #207

• Clustering of objects
  Divide the grid in subgrids and add an external additinal loop running over those subgrids (in parallel).
  But before executing this loop, all objects should be assigned to all subgrids in which they might be located.
  Thus, in each subgrid there are less objects to be checked reducing the calls of in.
  This might also be used in order to speed up the charge drift.

The URL of this package does not match that stored in METADATA.jl.
cc: @julia-tagbot[bot]

This should make proper use of the Unitful units as for example:

• understand numbers with units when plotting
• adapting the plot labels/ticks/units to the units given in the config file
• etc.

At the moment, we are simulating quite a number of detectors in our tests.
We should reduce that.

I guess besides the comparison to analytic solutions we only need to simulate 2 real detectors.
One for cartesian and one for cylindrical coordinates.

We should add then smaller, more specific tests for the geometry, e.g., primitives.
But we don't need a full simulation for each primitive.

Hello everyone,

When running the plot_electric_field(!) functions on a simulation with a cartesian geometry it does not show the expected behaviour. Instead, it outputs this
[ Info: 0 drifts are now being simulated...
┌ Warning: Attribute alias ylabel detected in the user recipe defined for the signature (::SolidStateDetectors.Plot_electric_field). To ensure expected behavior it is recommended to use the default attribute yguide.
└ @ Plots ~/.julia/packages/Plots/uCh2y/src/pipeline.jl:15
┌ Warning: Attribute alias xlabel detected in the user recipe defined for the signature (::SolidStateDetectors.Plot_electric_field). To ensure expected behavior it is recommended to use the default attribute xguide.
└ @ Plots ~/.julia/packages/Plots/uCh2y/src/pipeline.jl:15

And simply shows the electric potential with no field lines. The electric field calculation itself is fine to my knowledge.

Thanks!

Hello everyone,

I am wondering how the code actually deals with the impurity density profile and depletion. My current assumption is that the code assumes a standard pn junction to have a vanishingly small heavily doped side and assumes the 1D approximation where the depletion depth is independent of the heavily doped side. If this is the case, I understand that this is fine in Ge detectors because of the substantial crystal size. We would like to explicitly simulate both a p+ window (on n type bulk) and n+ implantation and pstop/pspray on the Ohmic side with varying thickness. This raises a number of questions:

• When defining the charge density model, I assume what is put in is instead the number density of the dopants. Is this correct? Does it assume everything is fully ionized? Does the sign matter in the config file if n/p type is specified?
• Is there an assumption of charge neutrality in the crystal when doing the calculations?
• Would it be correct to approximate these p+ and n+ layers as bulk semiconductor material with high 'charge' density in your program? How does the depletion work here? As an example, if I explicitly put a small layer of p+ Si on top of an n type crystal and apply a negative voltage to the p+ side, it seems like the p+ layer is undepleted when looking at the point types (yellow vs green), although the bulk is depleted. [This is also a more general question as to what each number of the point type stands for.]

Thanks!

Hi All,
I've constructed a 3 pixel planar detector and find that the field calculation is consuming 47GB of memory. I'm attaching my geometry file for reference. It's likely that I have something set up incorrectly. Any suggestions would be welcome.

simple_si_3hex.json.txt

We should change our definitions of all primitives to be centered around (0,0,0).
Like in GEANT4.

This also simplifies the rotation of objects (see #62).

It also comes in handy when we want to add some clustering of objects onto subregions of a grid
in order to parallelize and speed up the initialization of the boundary conditions on the grid for the SOR
which can take some time if there are to many objects and the grid becomes fine.
This clustering might also be used to speed up the charge drift in general (lesser lookups) (see #85 ).

• Using StructArrays (e.g. Path = Vector{CartesianPoint{T})
• Changing the naming of simulate! to something more specific
• Adjust simulate_waveforms for handling Tables
• Proper implementation of add_fano_noise
• Have a clear simulation chain: energy deposits -> charge deposits -> clustering + charge clouds -> drift -> waveforms
• Proper models for diffusion and self-repulsion and allowing the user to choose

I've been working on some more complex primitives that work within your CSG that might be of use for some users that have complex detector geometry that is hard to define implicitly. The detectors I've been simulating are for AGATA and as such are asymmetrical large-volume HPGe crystals with complex segmentation.

I struggled to create a good representation using your constructor so I wrote a primitive that can take any triangular mesh as an input (e.g. a standard stl, obj or ply file) and evaluate to the same level as your other primitives. The approach uses a generalised winding number for robust in-out segmentation. I've tested it with some simple geometry and it looks fine.

Here's a simple PoC using a modified version of your example, simple_mesh_simulation.pdf

Plotting works too however the 3D electrode plot takes some time to crunch through. Unfortunately for anything too complex the execution time and memory usage will crash the computer before completing, I've tried multithreading it but it's still too slow.

To solve this issue I decided to pre-evaluate the geometry on a voxelised grid and use that logic for geometry validation, the static nature of this method could limit spatial resolution however with small enough voxels I imagine true physical effects will be the limiting factor. Currently voxelisation is done externally in Python however I can likely implement a similar solution in Julia effectively using the code for the mesh primitive mentioned earlier.

Voxel evaluation is done on a regular grid to help speed up indexing, the approach works well however I'm not confident that the 'interesting' points are being properly passed. It's likely I can combine the two primitives to take in a standard mesh, voxelise it and determine the interesting points and then give that to the CSG.

Hopefully some of this is interesting to you.

It would be nice to allow for units in the user interface, especially for plotting, e.g.

plot(simulation.electric_field, xlims = (-4u"cm", 4u"cm"))

This also comes with nice unit handling in configuration files (see #182).

Hello,

We'd like to simulate a small part of the geometry of our Silicon detectors close to pixel edges, where we have substantial variation in the impurity density profile. Since the physical dimension of this region is very small compared to the total geometry, performing a regular simulation with extremely high grid density is not feasible. How would one best go about this in your code? I'm aware only of setting the boundary conditions to be set at infinity, is anything else included right now?

Thank you,
Leendert

Currently, only one single axis tick is inserted between two ticks if the electric potential difference, ΔΦ,
is above a certain limit, ΔΦ_max, independent of how large the difference is.
Thus, a lot of refinements steps are needed to get a good grid.

We should allow that multiple, N, ticks can be inserted between two ticks in one refinement step:
N = round(ΔΦ / ΔΦ_max)

By this we will also have to change the keywords regarding the grid and refinement in functions like
simulate! and calculate_electric_potential!.

Something like get_depletion_voltage(sim) that takes the Simulation, performs some kind of calculation and returns the depletion voltage in V.

Right the user is able to define the world in the config file.
Symmetries can be used there. For example mirror symmetries.

The issue is that the symmetries might differ for different things to simulate.
For example, the example coax detector has in general a periodicity of 60° in φ.
In addition, within this period, it has a mirror symmetry.
Thus, the electric potential only needs to be calculated on the interval [0°, 30°] (with reflecting boundary conditions in φ).

However, this is not possible for the weighting potentials of individual segments. Those, in this case, have
only a mirror symmetry but not a periodicity. Thus, they are calculated over 360°.
It could be reduced to [x°, x°+180°].

For this we need to add support for defining those symmetries in the configuration files for each individual contact.

ToDo:

• Support cylindrical coordinates

• Support Cartesian coordinates

• Support Depletion handling

• Add Documentation

• Add tests

• Test it on AMD GPUs

• How to choose a specific GPU device (-> using multiple GPUs in parallel)

• Add convergence check for GPU backend

Right now, the drift model assumes that the crystal axis <001> is parallel to the z-axis.

This issue is kinda related to #62.

Hello,

In our geometry we have a number of very thin layers compared to the overall dimensions of the system. One such thin layer is a pspray/pstop p+ layer which separates n+ pixels on n-type bulk. We would like to investigate what the behaviour of the system and pulse shapes is near pixel edges, which requires us to zoom into this region.

In doing preliminary testing, I tried to see what the code did when the physical size of a contact is much smaller than the anticipated grid spacing. As an example, I constructed a 10mm long tube with a 1um radius, with a grounded contact on one end and -1500V on the other, with a grid size of 30x30x30mm. Applying the initial conditions shows a large blot at the contact location with the correct voltage, which is fine. Performing the calculation with a number of refinements, however, a number of problems appear. Even with the convergence limit set to 1e-3, the process becomes extremely slow. More worryingly, it does not seem to reach convergence before stopping with no clear indication as to why. Finally, however, the potential obtained does not make sense, and the electric potential shows the tube having potentials of more than tens of thousands of positive voltage, with higher voltages obtained for larger grid spacings.

This raises a number of questions, some of which have been on our minds already a bit longer

• How is the convergence determined, and when does the algorithm decide to stop (and go to another refinement)? Is there a way to force it to continue?
• What does the 'value' stand for when the calculation is ongoing?
• Is this overshooting situation avoidable? Given that the initial state recovers the correct voltage for a very crude grid, I assume this is possible in principle
• Is it possible to have an adaptive (ie unequal spacings depending on the geometry in that location) grid?
• If not, is it possible for the user to define a custom grid with unequal spacings?

Thank you!
Leendert

We should implement a way to link other JSON files in one JSON config file to outsource larger fractions and make it easy to switch parts.

Status of the code: ready

To-Do:

• Function name decision:

scan_dict!(), deep_merge!(), recursive_merge!() ,... ?

• Folder in the source code
• Test for the function. Should there be a test added?
• Implementation in SSD. Suggestion in function:

parse_config_file(filename::AbstractString)::Dict{Any,Any}

at return

parsed_dict -> scan_dict!(parsed_dict)

• Documentation has to be written
• Should this be an internal or exported function?

I can do most of it, but I need answers to the points above from you.

I attached an example configuration with a jupyter notebook:

json-merger.zip

Hello,

When preparing toy examples for some of the previous issues posted today, I noticed that an asymmetry in the point types (and upon later inspection also on the electric potentials and every other quantity) was not a result of an asymmetric x-y geometry (like a hexagon), but also appears in the case of round tubes. As an example, I've attached the point types and electric potential of a simple tube in different planes.

big_tube.txt

Looking also at edges which are not parallel to x or y, there appear grid artefacts which can be tamed somewhat with increasing grid spacing, but remain visible.

In one case, this became particularly apparent

This is the central pixel of a hexagonal pixel ring. Every pixel contains roughly the same asymmetry, with a minimum along the x-axis through the center of the pixel and a butterfly structure around it.

Thank you!
Leendert

When I create a tube object with a translation and apply a cylindrical charge density with a radial gradient the charge density it is assumed the gradient starts from the z-axis. For example:

Is it possible to translate the gradient so that it is centered on the object instead?

The current implementation of Simulation{T} looks like this:

mutable struct Simulation{T <: SSDFloat} <: AbstractSimulation{T}
    detector::Union{SolidStateDetector{T}, Missing}
    q_eff_imp::Union{EffectiveChargeDensity{T}, Missing} # Effective charge coming from the impurites of the semiconductors
    q_eff_fix::Union{EffectiveChargeDensity{T}, Missing} # Fixed charge coming from fixed space charges, e.g. charged up surface layers
    ϵ_r::Union{DielectricDistribution{T}, Missing}
    point_types::Union{PointTypes{T}, Missing}
    electric_potential::Union{ElectricPotential{T}, Missing}
    weighting_potentials::Vector{Any}
    electric_field::Union{ElectricField{T}, Missing}

    charge_drift_model::Union{<:AbstractChargeDriftModel{T}, Missing}

    electron_drift_field::Union{ElectricField{T}, Missing}
    hole_drift_field::Union{ElectricField{T}, Missing}
end

with the default constructor being

function Simulation{T}() where {T <: SSDFloat}
    Simulation{T}(
        SolidStateDetector{T}(),
        missing,
        missing,
        missing,
        missing,
        missing,
        [missing],
        missing,
        ElectricFieldChargeDriftModel{T}(),
        missing,
        missing
    )
end

I would be in favor of replacing missing with something more type-stable. Why?
If simulation.electric_potential = missing, then calling plot(simulation.electric_potential) results in an error that missing cannot be converted to Float64 to be plotted rather than an error message that simulation.electric_potential is not yet calculated.

One way of doing this could be introducing new types, e.g. MissingPointTypes{T} <: PointTypes{T} and changing to point_types::PointTypes{T}. What would be your comments/thoughts on this?

I started refactoring our CSG. The core is finished. See branch refactorCSG.

This also concerns #62, #63, #84

The CSG is separated in a new submodule ConstructiveSolidGeometry inside SSD.

Rotation, translation and stretching/scaling are now implemented.

For now, I only added one primitive, the Tube.

I added a test script under /test/ConstructiveSolidGeometry/CSG_test.jl.
At the moment is more an example script.

We can start implementing our other primitives and functionality, like read-in, to it or adept to it.
But when adding something, please test for performance (type stability) like it is done in the test script.
A PR to this branch must have tests.

Guideline for Primitives

• Our primitives will be parametric types so that we can save computation time for special version of the primitives.
  E.g., in case of a full tube (cylinder) we don't have to check for rMin and the polar coordinate.

• Default constructors: All primitives should be constructed around (0,0,0). This makes rotations easier for the user. This is also the way how it is done in GEANT-4. E.g., Tube(radius = 1.0, height = 1.0) will create a Tube with an radius of 1 and going from -0.5 to +0.5 in z. Here, we will have to change (and break) our read-in our config files.

General remarks

• Branch refactorCSG: The old CSG is still inside the package. The new one is only included but not used yet.
  CartesianPoint and CylindricalPoint are not exported anymore to avoid conflicts as there exists now two version (for now). E.g.: SolidStateDetectors.CartesianPoint and SolidStateDetectors.ConstructiveSolidGeometry.CartesianPoint.

• The conversion of a CartesianPoint into a CylindricalPoint got changed.
  Previously: r = sqrt(pt.x * pt.x + pt.y * pt.y
  New: r = hypot(pt.x, pt.y) (special algorithm, see Julia documentation of hypot)

ToDo:

EDIT May 15th, 2021
The primitives have all been updated. Find a list of what still needs to be done at the end of this issue.

Clarify notation of the terms of impurity vs charge density and the naming of the corresponding functions and variables.

The field bulletize_PC in SigGen config files defines the point contact as spherical rather than as a tube.
This could be implemented using our Ellipsoid primitive.

Some of the functions are still left over from refactoring the JSON files.

One example is the bounding_box function, which returns the crystal_radius and crystal_length which is no longer stored in SolidStateDetector.

Also, we should look for functions that are defined but never called, e.g. from_internal_units.

These functions should be either removed or incorporated accordingly.

Based on the clarification of impurity and charge densities (#65), the redefinition of primitives (#84), the implementation of transformations (#62) and the refactoring of the CSG (#110), we need to update the config file format.

For this, we should decide on a long-term supported format, such that the user will have to update the config file format once (ideally with the help of conversion functions that SSD will provide).

When updating the configuration files to avoid deprecation warnings, I stumbled across a possible bug (?) in the painting algorithm.

Minimum working example

Configuration file: public_CGD_config_new.txt --> rename to public_CGD_config_new.yaml
(Sorry, but files that end with .yaml cannot be uploaded easily here)

using SolidStateDetectors
T = Float32

sim = Simulation{T}("public_CGD_config_new.yaml")
simulate!(sim, convergence_limit = 1e-6, refinement_limits = [0.2], verbose = false)
evt = Event([CartesianPoint{T}(0,2e-3,0)])
simulate!(evt, sim, Δt = 1e-9, max_nsteps = 10000)
signalsum = sum(broadcast(wf -> abs(wf.value[end]), evt.waveforms))

@info signalsum
@assert isapprox( signalsum, T(2), atol = 1e-3 )

It looks like the second contact does not get painted onto the grid.

Trouble shooting

The problem seems to be this geometry in the difference of the contact with id: 2:

- box:
    widths: [0.5, 10, 10]
    origin: 
      x: 4.75

This does not happen when widths: [0.6, 10, 10], when òrigin: 4.8, or when T = Float64.
So, I am assuming that this is connected to some kind of rounding error (?)

The geometry in the initial configuration file public_CGD_config.yaml also works perfectly and looks like this:

- box:
    x: 
      from: 4.5
      to: 5
    y: 10
    z: 10

Following #76, GR is now the standard plotting backend.
Some things need to be implemented accordingly:

1. Scalar potentials:

   • Equipotentials (fixed with #80)
   • Cross sections in z with nonuniform polar plots (fixed with #100)

2. Electric Field (#89):

   • Plot magnitude of the electric field
   • Plot electric field lines

3. Detectors:

   • Introduce consistent recipes for all primitives (Sphere could be plotted with PartialCircle if rotations are included)

Hello,

We are looking at the C-V behaviour of our geometries, but it is not entirely clear how to combine the information from point types and the capacitance calculation. Specifically, the capacitance continues decreasing even though (visually) all point types show up as green. As an example, taking a 2mm long tube with contacts on either end and varying the applied voltage and setting the impurity concentration to , you obtain the following C-V and d(1/C^2)/dV curves

From about 800V onwards, however, the point types show up as green, even though the slope of the C-V is non-zero at 1000V.

At 120V

At 780V

All electric potentials are calculated with depletion_handling=true on. In the naive analytical result one would expect the detector to be fully depleted around 600V, though I can understand edge effects playing a role here.

Thanks!
Leendert

big_tube.txt

Hello everyone,

It appears that currently the IO capabilities for HDF5 saving and loading are located in an unregistered package, LegendHDF5IO. Is it possible to integrate this directly into SSD, or if not bind the github repo's together?

Thanks!

The REQUIRE file could not be found.
cc: @oschulz

Hi,

I noticed a mistake in the hexagon primitive that was added recently and I thought it would be easier to fix through here than a full pull request.

In HexagonalPrism.jl line 76 should be
for φ in [deg2rad(30),deg2rad(90), deg2rad(150), deg2rad(210), deg2rad(270), deg2rad(330)]
Currently the plot recipe is rotated 30 degrees relative to the rest of the geometries functions.

For example,

versus

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.

Implement plot recipes for plotting the magnitude of the ElectricField and corresponding field lines.

Hello,

This might be more difficult to track down, but I've been experiencing segmentation faults over the past weeks when trying to plot quantities from SSD.

I have attached the full stack trace of the one that just occurred in the hope that you will be able to extract something useful.

julia> plot(s.electric_potential, x=0)
QWidget::paintEngine: Should no longer be called
QPainter::begin: Paint device returned engine == 0, type: 1
QWidget::paintEngine: Should no longer be called
QPainter::begin: Paint device returned engine == 0, type: 1

signal (11): Segmentation fault: 11
in expression starting at REPL[21]:0
_ZN13QMacCGContextC2EP8QPainter at /Users/leenderthayen/.julia/conda/3/lib/libQt5Gui.5.9.7.dylib (unknown line)
_ZN23QPixmapDropShadowFilter11qt_metacallEN11QMetaObject4CallEiPPv at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN8QToolBar10paintEventEP11QPaintEvent at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN11sipQToolBar10paintEventEP11QPaintEvent at /Users/leenderthayen/.julia/conda/3/lib/python3.7/site-packages/PyQt5/QtWidgets.so (unknown line)
_ZN7QWidget5eventEP6QEvent at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN8QToolBar5eventEP6QEvent at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN11sipQToolBar5eventEP6QEvent at /Users/leenderthayen/.julia/conda/3/lib/python3.7/site-packages/PyQt5/QtWidgets.so (unknown line)
_ZN19QApplicationPrivate13notify_helperEP7QObjectP6QEvent at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN12QApplication6notifyEP7QObjectP6QEvent at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN15sipQApplication6notifyEP7QObjectP6QEvent at /Users/leenderthayen/.julia/conda/3/lib/python3.7/site-packages/PyQt5/QtWidgets.so (unknown line)
_ZN16QCoreApplication15notifyInternal2EP7QObjectP6QEvent at /Users/leenderthayen/.julia/conda/3/lib/libQt5Core.5.9.7.dylib (unknown line)
_ZN14QWidgetPrivate14sendPaintEventERK7QRegion at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN14QWidgetPrivate10drawWidgetEP12QPaintDeviceRK7QRegionRK6QPointiP8QPainterP19QWidgetBackingStore at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN14QWidgetPrivate22paintSiblingsRecursiveEP12QPaintDeviceRK5QListIP7QObjectEiRK7QRegionRK6QPointiP8QPainterP19QWidgetBackingStore at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN14QWidgetPrivate22paintSiblingsRecursiveEP12QPaintDeviceRK5QListIP7QObjectEiRK7QRegionRK6QPointiP8QPainterP19QWidgetBackingStore at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN14QWidgetPrivate22paintSiblingsRecursiveEP12QPaintDeviceRK5QListIP7QObjectEiRK7QRegionRK6QPointiP8QPainterP19QWidgetBackingStore at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN14QWidgetPrivate10drawWidgetEP12QPaintDeviceRK7QRegionRK6QPointiP8QPainterP19QWidgetBackingStore at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN14QWidgetPrivate10scrollRectERK5QRectii at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN14QWidgetPrivate10scrollRectERK5QRectii at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN14QDesktopWidget11qt_metacallEN11QMetaObject4CallEiPPv at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN14QDesktopWidget11qt_metacallEN11QMetaObject4CallEiPPv at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN19QApplicationPrivate13notify_helperEP7QObjectP6QEvent at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN12QApplication6notifyEP7QObjectP6QEvent at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN15sipQApplication6notifyEP7QObjectP6QEvent at /Users/leenderthayen/.julia/conda/3/lib/python3.7/site-packages/PyQt5/QtWidgets.so (unknown line)
_ZN16QCoreApplication15notifyInternal2EP7QObjectP6QEvent at /Users/leenderthayen/.julia/conda/3/lib/libQt5Core.5.9.7.dylib (unknown line)
_ZN22QGuiApplicationPrivate18processExposeEventEPN29QWindowSystemInterfacePrivate11ExposeEventE at /Users/leenderthayen/.julia/conda/3/lib/libQt5Gui.5.9.7.dylib (unknown line)
_ZN22QWindowSystemInterface22sendWindowSystemEventsE6QFlagsIN10QEventLoop17ProcessEventsFlagEE at /Users/leenderthayen/.julia/conda/3/lib/libQt5Gui.5.9.7.dylib (unknown line)
_ZN22QWindowSystemInterface23flushWindowSystemEventsE6QFlagsIN10QEventLoop17ProcessEventsFlagEE at /Users/leenderthayen/.julia/conda/3/lib/libQt5Gui.5.9.7.dylib (unknown line)
qt_plugin_instance at /Users/leenderthayen/.julia/conda/3/plugins/platforms/libqcocoa.dylib (unknown line)
_ZN7QWindow10setVisibleEb at /Users/leenderthayen/.julia/conda/3/lib/libQt5Gui.5.9.7.dylib (unknown line)
_ZN14QWidgetPrivate8show_sysEv at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN14QWidgetPrivate11show_helperEv at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN7QWidget10setVisibleEb at /Users/leenderthayen/.julia/conda/3/lib/libQt5Widgets.5.9.7.dylib (unknown line)
_ZN14sipQMainWindow10setVisibleEb at /Users/leenderthayen/.julia/conda/3/lib/python3.7/site-packages/PyQt5/QtWidgets.so (unknown line)
_ZL17meth_QWidget_showP7_objectS0_ at /Users/leenderthayen/.julia/conda/3/lib/python3.7/site-packages/PyQt5/QtWidgets.so (unknown line)
_PyMethodDef_RawFastCallKeywords at /Users/leenderthayen/.julia/conda/3/lib/libpython3.7m.dylib (unknown line)
_PyCFunction_FastCallKeywords at /Users/leenderthayen/.julia/conda/3/lib/libpython3.7m.dylib (unknown line)
call_function at /Users/leenderthayen/.julia/conda/3/lib/libpython3.7m.dylib (unknown line)
_PyEval_EvalFrameDefault at /Users/leenderthayen/.julia/conda/3/lib/libpython3.7m.dylib (unknown line)
function_code_fastcall at /Users/leenderthayen/.julia/conda/3/lib/libpython3.7m.dylib (unknown line)
call_function at /Users/leenderthayen/.julia/conda/3/lib/libpython3.7m.dylib (unknown line)
_PyEval_EvalFrameDefault at /Users/leenderthayen/.julia/conda/3/lib/libpython3.7m.dylib (unknown line)
_PyEval_EvalCodeWithName at /Users/leenderthayen/.julia/conda/3/lib/libpython3.7m.dylib (unknown line)
_PyFunction_FastCallDict at /Users/leenderthayen/.julia/conda/3/lib/libpython3.7m.dylib (unknown line)
_PyObject_Call_Prepend at /Users/leenderthayen/.julia/conda/3/lib/libpython3.7m.dylib (unknown line)
method_call at /Users/leenderthayen/.julia/conda/3/lib/libpython3.7m.dylib (unknown line)
PyObject_Call at /Users/leenderthayen/.julia/conda/3/lib/libpython3.7m.dylib (unknown line)
macro expansion at /Users/leenderthayen/.julia/packages/PyCall/tqyST/src/exception.jl:95 [inlined]
#109 at /Users/leenderthayen/.julia/packages/PyCall/tqyST/src/pyfncall.jl:43 [inlined]
disable_sigint at ./c.jl:446 [inlined]
__pycall! at /Users/leenderthayen/.julia/packages/PyCall/tqyST/src/pyfncall.jl:42 [inlined]
_pycall! at /Users/leenderthayen/.julia/packages/PyCall/tqyST/src/pyfncall.jl:29
_pycall! at /Users/leenderthayen/.julia/packages/PyCall/tqyST/src/pyfncall.jl:11 [inlined]
#_#116 at /Users/leenderthayen/.julia/packages/PyCall/tqyST/src/pyfncall.jl:86 [inlined]
PyObject at /Users/leenderthayen/.julia/packages/PyCall/tqyST/src/pyfncall.jl:86
_display at /Users/leenderthayen/.julia/packages/Plots/uCh2y/src/backends/pyplot.jl:1376
display at /Users/leenderthayen/.julia/packages/Plots/uCh2y/src/output.jl:150
display at ./multimedia.jl:323
jl_apply at /Users/julia/buildbot/worker/package_macos64/build/src/./julia.h:1700 [inlined]
do_apply at /Users/julia/buildbot/worker/package_macos64/build/src/builtins.c:643
jl_f__apply at /Users/julia/buildbot/worker/package_macos64/build/src/builtins.c:657 [inlined]
jl_f__apply_latest at /Users/julia/buildbot/worker/package_macos64/build/src/builtins.c:693
#invokelatest#1 at ./essentials.jl:712 [inlined]
invokelatest at ./essentials.jl:711 [inlined]
print_response at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:161
print_response at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:146
unknown function (ip: 0x1d02551b5)
do_respond at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:729
unknown function (ip: 0x1d01f9413)
jl_apply at /Users/julia/buildbot/worker/package_macos64/build/src/./julia.h:1700 [inlined]
do_apply at /Users/julia/buildbot/worker/package_macos64/build/src/builtins.c:643
jl_f__apply at /Users/julia/buildbot/worker/package_macos64/build/src/builtins.c:657 [inlined]
jl_f__apply_latest at /Users/julia/buildbot/worker/package_macos64/build/src/builtins.c:693
#invokelatest#1 at ./essentials.jl:712 [inlined]
invokelatest at ./essentials.jl:711 [inlined]
run_interface at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/REPL/src/LineEdit.jl:2354
run_frontend at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:1055
run_repl at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:206
#764 at ./client.jl:383
jfptr_#764_8194.clone_1 at /Applications/Julia-1.4.app/Contents/Resources/julia/lib/julia/sys.dylib (unknown line)
jl_apply at /Users/julia/buildbot/worker/package_macos64/build/src/./julia.h:1700 [inlined]
do_apply at /Users/julia/buildbot/worker/package_macos64/build/src/builtins.c:643
jl_f__apply at /Users/julia/buildbot/worker/package_macos64/build/src/builtins.c:657 [inlined]
jl_f__apply_latest at /Users/julia/buildbot/worker/package_macos64/build/src/builtins.c:693
#invokelatest#1 at ./essentials.jl:712 [inlined]
invokelatest at ./essentials.jl:711 [inlined]
run_main_repl at ./client.jl:367
exec_options at ./client.jl:305
_start at ./client.jl:484
jfptr__start_2075.clone_1 at /Applications/Julia-1.4.app/Contents/Resources/julia/lib/julia/sys.dylib (unknown line)
true_main at /usr/local/bin/julia (unknown line)
main at /usr/local/bin/julia (unknown line)
Allocations: 17013236864 (Pool: 17013122412; Big: 114452); GC: 4127
Segmentation fault: 11

My working environment:

[a93c6f00] DataFrames v0.22.2
[8197267c] IntervalSets v0.5.2
[91a5bcdd] Plots v1.6.12
[92933f4c] ProgressMeter v1.4.1
[d330b81b] PyPlot v2.9.0
[bf2c0563] RadiationDetectorSignals v0.1.5
[71e43887] SolidStateDetectors v0.5.1 #master (https://github.com/JuliaPhysics/SolidStateDetectors.jl.git)
[276daf66] SpecialFunctions v1.2.1
[f3b207a7] StatsPlots v0.14.17
[1986cc42] Unitful v1.5.0
[69ce31c6] UpROOT v0.3.0

Find a list of things that need to be written/updated for the v0.6 documentation:

• How to define a ChargeDriftModel
• More information on the ADLChargeDriftModel (#144)
• How to define an ImpurityDensityModel (#140)
• How to create/structure config files (#136), especially concerning rotations (#62), crystal orientation (#161), and other stuff that got added
• How to split config files (#104)
• How to plot potentials and fields (ElectricPotential, ElectricField + lines, ...)
• Docstrings for exported functions
• Add a reference to the SSD publication (arXiv or JINST publication) to README.md (?)
  --> #175

Feel free to complete this list.

Seems f3b0488 was not enough to fix it 😢

I just found JuliaDocs/Documenter.jl#1277:
So we might also have to update our TagBot.yml, but it is not sure if it fixes the issue.

Hi,

It seems like when specifying a radial charge distribution gradient in a cartesian grid, it takes the 'r' argument and turns it into an x gradient instead.