Repository navigation
Fix gauge invariance in electric potential depletion-handling - #634
claudiaalvgar wants to merge 18 commits into
Conversation
|
Without having had a close look at the code, but looking at the description:
|
| # Seed the interior at minimum_applied_potential, not a hardcoded 0. | ||
| # A literal-0 start isn't gauge invariant -- it sits at a different | ||
| # point relative to the contacts depending on which one is | ||
| # grounded. Since depletion handling's clamp is bistable (issue | ||
| # #618), starting two gauge-equivalent solves from | ||
| # non-corresponding states can make them converge to different | ||
| # fixed points, even with a shift-invariant update rule. | ||
| potential = ismissing(potential_array) ? fill(minimum_applied_potential, size(grid)...) : potential_array |
There was a problem hiding this comment.
I'm not sure if this is what we want because minimum is also not invariant with respect to sign-flips.
There was a problem hiding this comment.
I fixed: the interior seed now uses gauge_ref_potential = mean(contact potentials), which is odd under sign flip (unlike min), and directly verified in tests. One thing to take into account is: :infinite's own boundary decay ended up staying on min, not mean because mean measurably corrupts the near-boundary field for realistic single-bias-direction detectors (confirmed on the isochrone-style geometry). So :infinite's truncation itself is not sign-flip-invariant but the cryostat test added shows that's an unavoidable property of any boundary approximation
| contact_bias_voltages::Vector{T} = if is_weighting_potential | ||
| T[0, 1] |
There was a problem hiding this comment.
Then, we might want to rename this to contact_potentials (instead of contact_bias_voltages) if we are also using this for weighting potentials.
| # Checks pn_junction_bit & update_bit, not bulk_bit. bulk_bit only | ||
| # gets set when ALL 26 neighbors of a point are also pn_junction -- | ||
| # so any point one grid step from a contact can never count as | ||
| # "bulk", even if it's genuinely undepleted. That's exactly where an | ||
| # undepleted sliver ends up right before full depletion: as bias | ||
| # approaches V_d, the undepleted region shrinks toward whichever | ||
| # contact it's pinching off against. So checking bulk_bit made | ||
| # is_depleted blind to tese cases. | ||
| !any(b -> (pn_junction_bit & b > 0) && (update_bit & b > 0) && (undepleted_bit & b > 0) && |
There was a problem hiding this comment.
Back in the days (#407), we included the bulk_bit check because we were getting individual grid points very close to the surface that would be wrongly marked as undepleted. How can we avoid those messing up the is_depleted result?
There was a problem hiding this comment.
In theory we could leave it as is as long as the grid size close to the surface is small enough. Sure it may lag the true depletion voltage by a bit, but would be more robust. Maybe we can quantify this effect with a reasonable grid size? I think for explorations of the Li profile and depletion voltage the grid size should be small enough by construction?
There was a problem hiding this comment.
I'm not sure we have to worry about the Li profile, because that part should be covered by the inactive_layer_bit part in the next line. I think this becomes a more relevant issue for coarse grids.
There was a problem hiding this comment.
I would also like to see some tests where all signs are swapped (bias voltage, impurity density) to cross-check if the electrical potential is also sign-flipped.
There was a problem hiding this comment.
I added new tests for this in the new commits
| using SolidStateDetectors | ||
| using Test |
There was a problem hiding this comment.
| using SolidStateDetectors | |
| using Test | |
| # This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT). | |
| using Test | |
| using SolidStateDetectors |
| apply_initial_state!(sim, ElectricPotential, Grid(sim)) | ||
| calculate_electric_potential!(sim, depletion_handling = false, refinement_limits = refinement_limits, use_nthreads = 4, verbose = false) |
There was a problem hiding this comment.
Read in test_utils.jl at the beginning and used the timed_ equivalent functions here (and elsewhere in this file).
| # undepleted region near the r=0 axis that used to be missed, which | ||
| # raises the true depletion voltage. So the old ~0.5% margin above | ||
| # `deplV` is no longer enough -- 15% above `deplV` is confirmed depleted. | ||
| sim.detector = SolidStateDetector(sim.detector, contact_id = id, contact_potential = ustrip(deplV * 1.15)) |
There was a problem hiding this comment.
I'm not sure i like this... In an ideal world, estimate_depletion_voltage should give a good estimate and not something that is 15% off..
There was a problem hiding this comment.
I think what's happening is that at deplV*1.005 (the old 0.5% margin), the detector is not depleted with the current fixes. This isn't a numerical issue in estimate_depletion_voltage, it's that deplV is an approximate linear estimate, while the true depletion threshold with depletion_handling on sits measurably higher once r0_handling correctly catches a previously-missed undepleted sliver near the axis
| @info signalsum | ||
| @test isapprox( signalsum, T(2), atol = 5e-3 ) | ||
| @test isapprox(timed_estimate_depletion_voltage(sim, verbose = false, check_for_depletion = false), T(-13.15)*u"V", atol = 1.0u"V") | ||
| # estimate_depletion_voltage's only(...) call assumes a single filtered |
There was a problem hiding this comment.
Why are we getting two candidates?
I can see this happening if there were one positive and one negative candidate, but for two negative candidates, we should always return the more extreme one..
There was a problem hiding this comment.
Should be fixed with #623, but I would still be interested why the changes in this PR result in two candidates instead of just one. Is this related to removing the bulk_bit check from is_depleted?
| # `dep_target` is matched analytically, via superposition, no solve, so | ||
| # no numerical noise. `dep_sim`, on the other hand, comes from a fresh SOR | ||
| # solve near the depletion threshold, which does carry the near-threshold | ||
| # stagnation noise described above. So the entire gap between the two | ||
| # numbers shows up in this comparison. | ||
| @test abs(dep_sim - dep_target) < 75u"V" | ||
| @test sim.detector.semiconductor.impurity_density_model != imp_model_before | ||
|
|
||
| # Re-run simulation in place and check depletion voltage matches again. This is a check that impurity_density_model and | ||
| # contact_potential where adapted correctly | ||
| # contact_potential where adapted correctly. | ||
| timed_calculate_electric_potential!(sim, refinement_limits = 0.01, depletion_handling = true) | ||
| @test abs(estimate_depletion_voltage(sim, check_for_depletion = false, verbose = false) - dep_sim) < 5u"V" | ||
| @test abs(estimate_depletion_voltage(sim, check_for_depletion = false, verbose = false) - dep_sim) < 75u"V" | ||
|
|
||
| # Finally, compare to fresh simulation which is changed manually | ||
| sim_fresh = Simulation{T}(joinpath(@__DIR__, "test_config_files/BEGe_01.yaml")) | ||
| sim_fresh.detector = SolidStateDetector(sim_fresh.detector, contact_id = id, contact_potential = bias_target) | ||
| sim_fresh.detector = SolidStateDetector(sim_fresh.detector, sim.detector.semiconductor.impurity_density_model) | ||
| timed_calculate_electric_potential!(sim_fresh, refinement_limits = 0.01, depletion_handling = true) | ||
| @test abs(estimate_depletion_voltage(sim_fresh, check_for_depletion = false, verbose = false) - dep_sim) < 5u"V" | ||
| @test abs(estimate_depletion_voltage(sim_fresh, check_for_depletion = false, verbose = false) - dep_sim) < 75u"V" |
There was a problem hiding this comment.
I'm also not sure I like the fact that the values are now 75V apart instead of 5V..
There was a problem hiding this comment.
I don't like this either, that's why I wanted to keep this as a draft. For Cartesian, unfixed code's gauges agree again above ~2030V too — not sure if that's a fluke or expected (agreement at high bias may just mean we've passed both gauges' thresholds, regardless of correctness — same pattern as r0-only Cylindrical above 2020V). I'll dig into the rest so we can reach a conclusion of whether this is a change we want to have or not
|
Thanks! |
| ϕmax - ϕmin < abs(U) ? Umin = U : Umax = U | ||
| end | ||
| end | ||
| bulk = findall((sim.point_types.data .& bulk_bit .> 0) .& (sim.point_types.data .& inactive_layer_bit .== 0)) |
There was a problem hiding this comment.
estimate_depletion_voltage still seems to be restricted to bulk_bits, whereas this restriction was removed in is_depleted.
Could this be why the results of estimate_depletion_voltage and manual is_depleted checks used to be only ~10V apart and are now 75V apart.
Can we adjust estimate_depletion_voltage to give more realistic estimates, e.g. by dealing differently with bulk_bits?
There was a problem hiding this comment.
I checked this and I think this is not the issue. BEGe_01 :infinite boundaries only had a 20mm buffer beyond the 40mm crystal. The Faraday-cage/sign-flip tests in test_gauge_invariance.jl show that :infinite's truncation-edge behaviour is only accurate once it has real separation from the region of interest. So the fix in my commits was to give BEGe_01.yaml the separation :infinite needs and all four comparisons now pass at 10V
| # pcs.potential lives in the V_ref = 0 (gauge-shifted) frame; shift back to real potentials once, here. | ||
| sim.electric_potential = ElectricPotential(ElectricPotentialArray(pcs) .+ pcs.gauge_ref_potential, grid) |
There was a problem hiding this comment.
Maybe this is something to add in ElectricPotentialArray instead of here?
There was a problem hiding this comment.
function ElectricPotentialArray(pcs::PotentialCalculationSetup{T, Cylindrical, 3, Array{T, 3}})::Array{T, 3} where {T}
pot::Array{T, 3} = Array{T, 3}(undef, size(pcs.grid))
for iz in axes(pot, 3)
irbz::Int = rbidx(iz)
for iφ in axes(pot, 2)
irbφ::Int = iφ + 1
idxsum::Int = iz + iφ
for ir in axes(pot, 1)
irbr::Int = ir + 1
rbi::Int = iseven(idxsum + ir) ? rb_even::Int : rb_odd::Int
pot[ir, iφ, iz] = pcs.potential[ irbz, irbφ, irbr, rbi ]
end
end
end
for iz in axes(pot, 3)
p_r0 = mean(pot[1,:,iz])
pot[1,:,iz] .= p_r0
end
# pcs.potential lives in the V_ref = 0 (gauge-shifted) frame; shift back to real potentials once, here.
return pot .+ pcs.gauge_ref_potential
end
86316d2 to
493883f
Compare
| # `deplV` is no longer enough -- 15% above `deplV` is confirmed depleted. | ||
| sim.detector = SolidStateDetector(sim.detector, contact_id = id, contact_potential = ustrip(deplV * 1.15)) |
There was a problem hiding this comment.
Can we somehow address to make estimate_depletion_voltage better in this sense? 15% tolerance is pretty bad in my opinion..
There was a problem hiding this comment.
@claudiaalvgar could you check if with the updates to estimate_depletion_voltage in PR #639, the two methods agree again (and we don't have to increase the value from estimate_depletion_voltage by 15% to make this work).
| @test (timed_estimate_depletion_voltage(sim, verbose = false, check_for_depletion = false); true) | ||
| id = SolidStateDetectors.determine_bias_voltage_contact_id(sim.detector) | ||
| sim.detector = SolidStateDetector(sim.detector, contact_id = id, contact_potential = T(-5.0)) | ||
| timed_calculate_electric_potential!(sim, depletion_handling = true) | ||
| @test !is_depleted(sim.point_types) | ||
|
|
||
| sim.detector = SolidStateDetector(sim.detector, contact_id = id, contact_potential = T(-11.0)) | ||
| timed_calculate_electric_potential!(sim, depletion_handling = true) | ||
| @test is_depleted(sim.point_types) |
There was a problem hiding this comment.
Also, what used to be a depletion voltage of -13.15V is now between -5V and -11V. Do we understand this, and does this make sense?
There was a problem hiding this comment.
-13.15V came from the previous code, where :infinite always decayed toward a hardcoded 0 on both the real solve and the weighting-potential solve used in estimate_depletion_voltage so no mismatch. This PR's gauge-invariance change makes the real solve's :infinite decay target depend on the bias sign, while the weighting potential's stays fixed at 0. For this detector's negative bias, that mismatch corrupts estimate_depletion_voltage to ≈-70 to -75V (in both Float32 and Float64). :reflecting's boundary formula is a plain mirror-copy of the neighbouring grid cell with no reference-potential term at all, so running this same code with :reflecting instead gives ≈-12.16V/-12.19V, close to the old -13.15V
| @@ -17,8 +17,8 @@ grid: | |||
| left: periodic | |||
| right: periodic | |||
| z: | |||
| from: -20 | |||
| to: 60 | |||
| from: -300 | |||
| to: 300 | |||
There was a problem hiding this comment.
I know that these changes were done to move the :infinite boundary further away from the crystal. But is there any way that we can also get a reasonable boundary condition without having to extend the world this much?
There was a problem hiding this comment.
The enlargement here is compensating for the same weighting-potential mismatch inside adjust_impurity_and_electric_potential_to_match_depletion!/adjust_bias_and_electric_potential!, which pushes dep_sim ≈30V off target on that small world. For now, the validated fix could be using :reflecting : no reference-potential term, so it can't develop this mismatch
|
I ran the combined PR634+639 (julia version 1.11.5): on the real-detector config, Float32 now gives 2030.11 V against Float64's 2030.006 V — so PR639's fix holds fully when combined with PR634, and the old 15% safety margin isn't needed anymore; on the BEGe config, though, the U_est/U_alt comparison in test_depletion.jl is unmoved by adding PR639 — both still fail the <5V check, but this failure would also be true if we don't consider PR634 nor PR639 and use Float64 instead of Float32. I put the results on the table attached for future reference |
029c31f to
3af9ea4
Compare
|
Update: I rebased this branch since the changes were needed for these fixes to apply cleanly. PR634 added |
…nferring it from the contacts
|
Update: explicit potential of the surroundings for The rule from my last update (
Changes
Results (InvertedCoax, orig/swap mismatches in bold)
Cylindrical (boundary on r_max, z_min, z_max)
Cartesian (InvertedCoax on an 80 × 80 × 100 mm box, boundary on all six faces)
|
…y works correctly
This branch has not been deployed
3 participants
Some issues in the depletion-handling / SOR machinery all shared the same root cause: an implicit, unshiftable "0 V" reference baked into several places that should instead reference the detector's actual applied potentials. Fixed together because they're only fully verifiable together, the r0-axis fix alone still leaves gauge-dependent mismatches near threshold (see table below).
BoundaryConditionsCartesian.jlThe
:infiniteboundary condition sets edge cells to decay toward literal 0, which is correct only when the lowest applied potential (minimum_applied_potential) equals 0V. Under gauge shifts, this decayed toward wrong values. Since the potential array itself now lives in thegauge_ref_potential-shifted frame (established once at setup, not per-iteration), reconciling it with:infinite's actual target requires a second shift:Δ = minimum_applied_potential - gauge_ref_potential, applied to just the edge margin before the formula runs and undone after:n → n - Δ → f·(n - Δ) → f·(n - Δ) + Δ. This runs every SOR iteration, on narrow edge margins only.BoundaryConditionsCylindrical.jlSame fix, applied to the r-axis and z-axis dispatchers for cylindrical grids, for the same reason.
PotentialCalculationSetupCartesian3D.jl/PotentialCalculationSetupCylindrical.jlTwo independent fixes each:
zeros(T, ...)tofill(gauge_ref_potential, ...)— a new struct field, the mean of the applied contact potentials (0for weighting-potential solves). An earlier version of this fix usedfill(minimum_applied_potential, ...)instead, butminisn't antisymmetric under a simultaneous sign flip of bias and impurity density (minandmaxswap contacts), so two sign-mirrored, gauge-equivalent configurations could still seed inconsistently and converge to different fixed points.gauge_ref_potential(mean) is exactly antisymmetric under that flip, fixing this.minimum_applied_potentialis retained separately as:infinite's own decay target (see above) — the two references now serve different purposes and are not the same value.is_weighting_potentialguard oncontact_potentials(renamed fromcontact_bias_voltages): Weighting-potential solves use the 0/1 convention (target = 1, rest = 0), not the detector's real bias voltages. Without this guard,minimum_applied_potential/gauge_ref_potentialwould derive from real (possibly large) bias voltages even during weighting-potential solves, seeding incorrect values.Also added:
estimate_depletion_voltage(Depletion.jl) now combines all filtered depletion-voltage candidates viamax/mininstead of assuming exactly one viaonly(...), fixing a crash on detectors with multiple candidates. The "shift back to real units" step for the electric potential now lives insideElectricPotentialArray(both coordinate systems) rather than duplicated at eachSimulation.jlcall site.CPU_innerloop.jlr0_handling_depletion_handling(Cylindrical) read an unwritten "ghost" neighbour at the r=0 axis that silently defaults to 0, acting as an implicit depletion floor regardless of gauge. Now duplicates the r-right (radially outward) neighbour into that slot instead — a value that moves with the actual solution rather than a fixed external constant.convergence.jlPurely diagnostic, no solved values change: labels why the SOR loop exited (
:converged/:stagnated/:max_iterations) and warns once (maxlog=1, separate_ids per cause) when stagnation is due to near-threshold depletion bistability vs. hitting floating-point precision. Near-threshold points can genuinely flip classification between iterations, so tight-tolerance convergence is not always reachable, by design.PointTypes.jlis_depletedrequiredbulk_bit(all 26 neighbours alsopn_junction_bit), which structurally excludes any point one grid-step from a contact — exactly where a shrinking undepleted sliver ends up right before true full depletion. Switched topn_junction_bit & update_bit, with no such neighbour requirement — the same pairget_active_volumealready uses elsewhere in this file for the equivalent question.Tests
test_gauge_invariance.jl(new): Validates fixes against closed-form analytical depletion voltages across ten testsets covering planar slabs, annular coax, and cylindrical geometries with various boundary conditions — including direct pointwise field-comparison checks (epotΔ ≈ epot0 + Δ) under bothdepletion_handlingvalues, which probe the mean-seed fix directly rather than only throughis_depleted.test_depletion.jl: Updatedr0_handling_depletion_handlingunit test to corrected output. Tolerances were initially loosened to account for observed near-threshold noise, but that noise was subsequently root-caused toBEGe_01.yaml's:infiniteboundaries sitting too close to the crystal (a 20mm buffer around a 40mm crystal).test_real_detectors.jl: Inverted Coax's post-depletion margin widened (0.5% → 15%) since the r0-axis fix now correctly detects a small undepleted region near the axis that used to be missed, raising the true depletion voltage. Hexagon's depletion-voltage check replaced with an is_depleted-based bracket: confirmed undepleted at -5.0V and depleted at -11.0V.test_refine_existing_potential.jl: grid-size assertions updated.Gauge-invariance sweep, Cylindrical (InvertedCoax)
T = Float64 refinement = [0.2, 0.1, 0.05]
:infinite boundary condition
All-fixes matches ground truth at every single voltage tested, including exactly at and above the true threshold (2015.76V–2018V), where r0-only alone still leaves 3 gauge-dependent mismatches — direct evidence that the
:infinite-boundary and interior-seed fixes are independently necessary, not redundant with the r0-axis fix.Gauge-invariance sweep, Cartesian (InvertedCoax)
No r0-only column here: Cartesian grids never dispatch through
r0_handling_depletion_handlingat all, so that fix is structurally irrelevant to this case — only the:infinite-boundary and interior-seed fixes can do anything for Cartesian.T = Float64 refinement = [0.2, 0.1, 0.05]
:infinite boundary condition
All-fixes matches at every point. Confirms the Cylindrical result generalises: on a coordinate system where the r0 fix contributes nothing at all, the remaining two fixes alone are sufficient to fully restore gauge invariance here.
Holding the algorithm fixed and varying only the boundary type shows the near-threshold transition point itself moving by several volts:
:reflectingtransitions around 2022-2024V,:fixedaround 2024-2025V, and:infinite(table above) around 2028-2029V.Open questions / still needs verification
For detectors using only
:infiniteboundaries,estimate_depletion_voltage's numeric answer is sensitive to how much separation the grid gives the boundary from the crystal. Confirmed directly: enlargingBEGe_01.yaml's grid (20mm→300mm buffer) reduced a comparable instability from ~50V down to under 1V. This is still an open, structural limitation of:infinite(see the Faraday-cage test intest_gauge_invariance.jl) — not something a code fix removes, since the correct open-boundary reference is geometry/capacitance-dependent, not derivable from contact potentials alone. Giving:infiniteenough separation mitigates it; nothing eliminates it outright.