Add adapt_to_depletion_and_bias_voltage!

src/Simulation/Depletion.jl

@@ -204,6 +204,76 @@ function _find_depletion_voltage_candidates(ϕᵨ::AbstractArray{T, 3}, ϕᵥ::A
     minimum(Umin), maximum(Umax)
 end
 
+"""
+    adapt_to_depletion_and_bias_voltage!(sim::Simulation{T}, dep::RealQuantity, bias::RealQuantity;
+        contact_id::Int = determine_bias_voltage_contact_id(sim.detector),
+        verbose::Bool = true,
+        kwargs...) where {T <: AbstractFloat}
+
+Adapts a [`Simulation`](@ref) in place so that it matches a target depletion voltage `dep`
+and bias voltage `bias`, **without re-solving the field**.
+
+This works by exploiting the linearity of the electric potential in both the impurity density
+and the applied bias. The impurity density model is rescaled by `f = dep / dep_sim`, where `dep_sim`
+is the current depletion voltage estimated via [`estimate_depletion_voltage`](@ref), and the bias
+contact potential is swapped to `bias`. The stored `sim.electric_potential` is updated accordingly via
+superposition with the [`WeightingPotential`](@ref) of the bias contact, so no new field calculation
+is performed.
+
+`dep` and `bias` must share the same (non-zero) sign, and `bias` must exceed `dep` in magnitude
+(the detector must be over-depleted at the operating voltage).
+
+## Arguments
+* `sim::Simulation{T}`: [`Simulation`](@ref) to be adapted in place.
+* `dep::RealQuantity`: Target depletion voltage. If no units are given, this value is parsed in units of `$(internal_voltage_unit)`.
+* `bias::RealQuantity`: Target operating (bias) voltage. If no units are given, this value is parsed in units of `$(internal_voltage_unit)`.
+
+## Keywords
+* `contact_id::Int`: The `id` of the [`Contact`](@ref) at which the bias voltage is applied.
+    The default is determined automatically via `determine_bias_voltage_contact_id(sim.detector)`.
+* `verbose::Bool = true`: Activate or deactivate additional info output. Default is `true`.
+
+Additional `kwargs...` are passed on to [`estimate_depletion_voltage`](@ref).
+
+## Example
+```julia
+using SolidStateDetectors
+sim = Simulation(SSD_examples[:InvertedCoax])
+calculate_electric_potential!(sim)
+adapt_to_depletion_and_bias_voltage!(sim, 1000u"V", 1500u"V")
+```
+
+!!! note
+    The accuracy of the result depends on the precision of the initial simulation, since the
+    impurity density is rescaled rather than re-solved.
+
+See also [`estimate_depletion_voltage`](@ref).
+"""
+function adapt_to_depletion_and_bias_voltage!(sim::Simulation{T}, dep::RealQuantity, bias::RealQuantity;
+    contact_id::Int = determine_bias_voltage_contact_id(sim.detector),
+    verbose::Bool = true,
+    kwargs...) where {T <: AbstractFloat}
+    if ismissing(sim.weighting_potentials[contact_id]) || sim.weighting_potentials[contact_id].grid != sim.electric_potential.grid
+        _adapt_weighting_potential_to_electric_potential_grid!(sim, contact_id)
+    end
+    dep::T = _parse_value(T, dep, internal_voltage_unit)
+    bias::T = _parse_value(T, bias, internal_voltage_unit)
+    @assert dep * bias > 0 "The depletion voltage ($(dep)$(internal_voltage_unit)) and operating voltage ($(bias)$(internal_voltage_unit)) must have the same (non-zero) sign."
+    @assert abs(bias) > abs(dep) "The operating voltage ($(bias)$(internal_voltage_unit)) must exceed the depletion voltage ($(dep)$(internal_voltage_unit)) in magnitude (the detector must be over-depleted)."
+    dep_sim = _parse_value(T, estimate_depletion_voltage(sim; contact_id = contact_id, verbose = verbose, kwargs...), internal_voltage_unit)
+    f = T(dep/dep_sim)
+    if verbose
+        @info """Adapting `sim.electric_potential`, 
+        scaling `impurity_density_model` by $f, 
+        and swapping in the contact potential on `sim.detector.contacts[$contact_id]` 
+        such that the simulation matches the given depletion and bias voltages."""
+    end
+    ϕV = sim.weighting_potentials[contact_id].data   
+    sim.electric_potential.data .= f * sim.electric_potential.data .+ (bias - f * sim.detector.contacts[contact_id].potential) .* ϕV
+    sim.detector = SolidStateDetector(sim.detector, contact_id = contact_id, contact_potential = bias)
+    sim.detector = SolidStateDetector(sim.detector, f*sim.detector.semiconductor.impurity_density_model)
+    nothing
+end
 #=
 """
     old_estimate_depletion_voltage( sim::Simulation{T}, contact_id::Int, field_sim_settings = (verbose = true,))::T

src/SolidStateDetectors.jl

@@ -76,7 +76,7 @@ export ElectricFieldChargeDriftModel, ADLChargeDriftModel, ADL2016ChargeDriftMod
 export ConstantImpurityDensity, LinearImpurityDensity, ThermalDiffusionLithiumDensity, PtypePNJunctionImpurityDensity
 export LinBouleImpurityDensity, ParBouleImpurityDensity, LinExpBouleImpurityDensity, ParExpBouleImpurityDensity, SplineBouleImpurityDensity
 export NoChargeTrappingModel, BoggsChargeTrappingModel, ConstantLifetimeChargeTrappingModel, CombinedChargeTrappingModel
-export get_active_volume, is_depleted, estimate_depletion_voltage, set_point_type_depletion_handling!
+export get_active_volume, is_depleted, estimate_depletion_voltage, set_point_type_depletion_handling!, adapt_to_depletion_and_bias_voltage!
 export calculate_stored_energy, calculate_mutual_capacitance, calculate_capacitance_matrix
 export simulate_waveforms
 export run_geant4_simulation

test/test_depletion.jl

@@ -32,4 +32,29 @@ T = Float32
     @test_throws Exception estimate_depletion_voltage(sim, -10, 0, tolerance = 20)
     @test_throws Exception estimate_depletion_voltage(sim, 0u"kg", 20u"kg")
     @test_logs (:info,) (:info,) (:warn, r".*not in the specified range.*") estimate_depletion_voltage(sim, U_est/3, 0)
+
+    # `adapt_to_depletion_and_bias_voltage!` rescales the impurity density and swaps in a new
+    # contact potential so the simulation matches a target depletion voltage `dep` and
+    # bias voltage `bias`, without re-solving the field. Check the round-trip: after adapting,
+    # the estimated depletion voltage should be ≈ `dep` and the bias contact should sit at `bias`.
+    dep_target = 2000u"V"
+    bias_target  = 2500u"V"
+    imp_model_before = sim.detector.semiconductor.impurity_density_model
+    adapt_to_depletion_and_bias_voltage!(sim, dep_target, bias_target, check_for_depletion = false, verbose = false)
+    @test sim.detector.contacts[id].potential == SolidStateDetectors._parse_value(T, bias_target, SolidStateDetectors.internal_voltage_unit)
+    dep_sim = estimate_depletion_voltage(sim, check_for_depletion = false, verbose = false)
+    @test abs(dep_sim - dep_target) < 5u"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
+    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"
+
+    # 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"
 end