JuliaPhysics / JuliaPhysics/SolidStateDetectors.jl
Electric-field gradient becomes first-order accurate at abrupt grid-spacing changes
- Dominant language
- Julia
- Stars
- 178
- Forks
- 59
- Avg merge
- 1d 20h
- Merged PRs (30d)
- 4
Description
(Found in an independent external review of PR #623 — commissioned by @oschulz, performed by an AI reviewer — and verified against the code. Pre-existing behavior, not introduced by that PR.)
`get_electric_field_from_potential` (`src/ElectricField/ElectricField.jl`, both the cylindrical and Cartesian versions, all three axes) computes interior derivatives as the *unweighted* average of the forward and backward difference quotients:
```julia
er = ( Δp_1/d_1 + Δp_2/d_2 ) / 2
```
This is second-order accurate only for equal spacing `d_1 == d_2`. On unequal spacing the correct three-point (quadratic-interpolation) derivative weights each one-sided slope by the *opposite* interval:
```
f'(x_i) ≈ d_2/(d_1+d_2) * Δp_1/d_1 + d_1/(d_1+d_2) * Δp_2/d_2
```
Counterexample: `V = x²` at `x_i = 0` with spacings `d_2 = h₋`, `d_1 = h₊`. The exact derivative is 0; the implemented average gives `(h₊ − h₋)/2`. In general the local error is `f''(x_i)/4 · (h₊ − h₋)` — first order in the spacing jump.
Why it matters here: SSD's grids are adaptively refined, so spacing jumps (up to `max_distance_ratio = 5` by default) sit exactly where refinement put them — near contacts and in strong-field regions. The electric field, and hence drift velocities and pulse shapes, are biased around those transitions. The effect is bounded by the curvature times the spacing jump, so smooth regions are fine, but refinement boundaries cut through the most physics-relevant zones.
Suggested fix and validation:
1. Replace the plain average by the weighted form above (identical results on uniform spacing, so reference comparisons on uniform grids stay bit-compatible).
2. Check the interplay with the boundary-point special cases (the `e *= 2` one-sided corrections next to non-updated points) — those are a separate heuristic and should be revisited with the same stencil logic.
3. Add manufactured-solution tests: `V = a + b·x` must be exact on any grid; `V = x²` on a deliberately nonuniform axis must give `E = -2x` at interior ticks (this discriminates the current implementation); for cylindrical, `V = r^m cos(mφ)` on refined `r`/`φ` axes.
Changing this alters field values on every adaptively refined simulation (at the refinement transitions), so it should be validated against the example-detector references and probably released with a note.
Contributor guide
No contributing guide indexed for this repository
Research direction
Start in src/ElectricField/ElectricField.jl at get_electric_field_from_potential, covering both cylindrical and Cartesian implementations and all three axes. Review the existing e *= 2 boundary-point cases, then add manufactured-solution coverage for linear and quadratic potentials on nonuniform grids, plus the cylindrical case. Done means weighted interior derivatives, validated example-detector references, and an appropriate release note.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- julia
- Domain
- backend
- Issue type
- Bug
- Difficulty
- 4/5
- Estimated time
- 3-5 days
- Activity status
- Quiet
- Clarity
- Mostly clear
- Newbie friendliness
- 52/100