Skip to content

Commit 6eec967

Browse files
lmh91lmh91
authored
Merge pull request #206 from JuliaPhysics/mt_initial_state
Improve application of boundary conditions
2 parents 16264e8 + 1cf0dba commit 6eec967

File tree

5 files changed

+181
-170
lines changed

5 files changed

+181
-170
lines changed

‎src/PotentialSimulation/PotentialSimulationSetups/PotentialSimulationSetupRBCartesian3D.jl

+74-65
Original file line numberDiff line numberDiff line change
@@ -1,59 +1,53 @@
1-
function set_point_types_and_fixed_potentials!(point_types::Array{PointType, N}, potential::Array{T, N},
2-
grid::Grid{T, N, Cartesian}, det::SolidStateDetector{T}; weighting_potential_contact_id::Union{Missing, Int} = missing,
3-
not_only_paint_contacts::Val{NotOnlyPaintContacts} = Val{true}(),
4-
paint_contacts::Val{PaintContacts} = Val{true}())::Nothing where {T <: SSDFloat, N, NotOnlyPaintContacts, PaintContacts}
1+
function set_passive_or_contact_points(point_types::Array{PointType, 3}, potential::Array{T, 3},
2+
grid::CartesianGrid3D{T}, obj, pot::T, use_nthreads::Int = 1) where {T}
3+
if !isnan(pot)
4+
@onthreads 1:use_nthreads for iz in workpart(axes(potential, 3), 1:use_nthreads, Base.Threads.threadid())
5+
@inbounds for iy in axes(potential, 2)
6+
for ix in axes(potential, 1)
7+
pt::CartesianPoint{T} = CartesianPoint{T}( grid.axes[1].ticks[ix], grid.axes[2].ticks[iy], grid.axes[3].ticks[iz] )
8+
if pt in obj
9+
potential[ ix, iy, iz ] = pot
10+
point_types[ ix, iy, iz ] = zero(PointType)
11+
end
12+
end
13+
end
14+
end
15+
end
16+
nothing
17+
end
518

6-
axx::Vector{T} = grid.axes[1].ticks
7-
axy::Vector{T} = grid.axes[2].ticks
8-
axz::Vector{T} = grid.axes[3].ticks
19+
function set_point_types_and_fixed_potentials!(point_types::Array{PointType, 3}, potential::Array{T, 3},
20+
grid::CartesianGrid3D{T}, det::SolidStateDetector{T};
21+
weighting_potential_contact_id::Union{Missing, Int} = missing,
22+
use_nthreads::Int = Base.Threads.nthreads(),
23+
not_only_paint_contacts::Val{NotOnlyPaintContacts} = Val{true}(),
24+
paint_contacts::Val{PaintContacts} = Val{true}())::Nothing where {T <: SSDFloat, NotOnlyPaintContacts, PaintContacts}
925

10-
for iz in axes(potential, 3)
11-
z::T = axz[iz]
12-
for iy in axes(potential, 2)
13-
y::T = axy[iy]
26+
@onthreads 1:use_nthreads for iz in workpart(axes(potential, 3), 1:use_nthreads, Base.Threads.threadid())
27+
@inbounds for iy in axes(potential, 2)
1428
for ix in axes(potential, 1)
15-
x::T = axx[ix]
16-
pt::CartesianPoint{T} = CartesianPoint{T}( x, y, z )
17-
18-
if !ismissing(det.passives)
19-
for passive in det.passives
20-
if !isnan(passive.potential)
21-
if pt in passive
22-
potential[ ix, iy, iz ] = if ismissing(weighting_potential_contact_id)
23-
passive.potential
24-
else
25-
0
26-
end
27-
point_types[ ix, iy, iz ] = zero(PointType)
28-
end
29-
end
30-
end
31-
end
29+
pt::CartesianPoint{T} = CartesianPoint{T}( grid.axes[1].ticks[ix], grid.axes[2].ticks[iy], grid.axes[3].ticks[iz] )
3230
if in(pt, det.semiconductor)
3331
point_types[ ix, iy, iz ] += pn_junction_bit
3432
end
35-
if NotOnlyPaintContacts
36-
for contact in det.contacts
37-
if pt in contact
38-
potential[ ix, iy, iz ] = if ismissing(weighting_potential_contact_id)
39-
contact.potential
40-
else
41-
contact.id == weighting_potential_contact_id ? 1 : 0
42-
end
43-
point_types[ ix, iy, iz ] = zero(PointType)
44-
end
45-
end
46-
end
4733
end
4834
end
35+
end
36+
if !ismissing(det.passives)
37+
for passive in det.passives
38+
pot::T = ismissing(weighting_potential_contact_id) ? passive.potential : zero(T)
39+
set_passive_or_contact_points(point_types, potential, grid, passive.geometry, pot, use_nthreads)
40+
end
41+
end
42+
if NotOnlyPaintContacts
43+
for contact in det.contacts
44+
pot::T = ismissing(weighting_potential_contact_id) ? contact.potential : contact.id == weighting_potential_contact_id
45+
set_passive_or_contact_points(point_types, potential, grid, contact.geometry, pot, use_nthreads)
46+
end
4947
end
5048
if PaintContacts
5149
for contact in det.contacts
52-
pot::T = if ismissing(weighting_potential_contact_id)
53-
contact.potential
54-
else
55-
contact.id == weighting_potential_contact_id ? 1 : 0
56-
end
50+
pot::T = ismissing(weighting_potential_contact_id) ? contact.potential : contact.id == weighting_potential_contact_id
5751
fs = ConstructiveSolidGeometry.surfaces(contact.geometry)
5852
for face in fs
5953
paint!(point_types, potential, face, contact.geometry, pot, grid)
@@ -63,10 +57,29 @@ function set_point_types_and_fixed_potentials!(point_types::Array{PointType, N},
6357
nothing
6458
end
6559

60+
function fill_ρ_and_ϵ!(ϵ::Array{T}, ρ_tmp::Array{T}, q_eff_fix_tmp::Array{T},
61+
::Type{Cartesian}, mpz::Vector{T}, mpy::Vector{T}, mpx::Vector{T}, use_nthreads::Int, obj) where {T}
62+
@inbounds begin
63+
@onthreads 1:use_nthreads for iz in workpart(axes(ϵ, 3), 1:use_nthreads, Base.Threads.threadid())
64+
for iy in axes(ϵ, 2)
65+
for ix in axes(ϵ, 1)
66+
pt::CartesianPoint{T} = CartesianPoint{T}(mpx[ix], mpy[iy], mpz[iz])
67+
if pt in obj
68+
ρ_tmp[ix, iy, iz]::T, ϵ[ix, iy, iz]::T, q_eff_fix_tmp[ix, iy, iz]::T = get_ρ_and_ϵ(pt, obj)
69+
end
70+
end
71+
end
72+
end
73+
end
74+
nothing
75+
end
6676

67-
function PotentialSimulationSetupRB(det::SolidStateDetector{T}, grid::CartesianGrid3D{T}, medium::NamedTuple = material_properties[materials["vacuum"]],
68-
potential_array::Union{Missing, Array{T, 3}} = missing; weighting_potential_contact_id::Union{Missing, Int} = missing,
69-
sor_consts::T = T(1), not_only_paint_contacts::Bool = true, paint_contacts::Bool = true ) where {T}#::PotentialSimulationSetupRB{T} where {T}
77+
function PotentialSimulationSetupRB(det::SolidStateDetector{T}, grid::CartesianGrid3D{T},
78+
medium::NamedTuple = material_properties[materials["vacuum"]],
79+
potential_array::Union{Missing, Array{T, 3}} = missing;
80+
weighting_potential_contact_id::Union{Missing, Int} = missing,
81+
use_nthreads::Int = Base.Threads.nthreads(),
82+
sor_consts::T = T(1), not_only_paint_contacts::Bool = true, paint_contacts::Bool = true ) where {T}
7083
is_weighting_potential::Bool = !ismissing(weighting_potential_contact_id)
7184

7285
@inbounds begin
@@ -151,9 +164,6 @@ function PotentialSimulationSetupRB(det::SolidStateDetector{T}, grid::CartesianG
151164
(grid_boundary_factor_z_left, grid_boundary_factor_z_right))
152165
end
153166

154-
# detector_material_ϵ_r::T = det.material_detector.ϵ_r
155-
# environment_material_ϵ_r::T = det.material_environment.ϵ_r
156-
157167
contact_bias_voltages::Vector{T} = if length(det.contacts) > 0
158168
T[contact.potential for contact in det.contacts]
159169
else
@@ -165,18 +175,14 @@ function PotentialSimulationSetupRB(det::SolidStateDetector{T}, grid::CartesianG
165175
depletion_handling_potential_limit::T = -bias_voltage
166176
sor_consts = [sor_consts]
167177

168-
ϵ = Array{T, 3}(undef, length(mpx), length(mpy), length(mpz))
169-
ρ_tmp = Array{T, 3}(undef, length(mpx), length(mpy), length(mpz))
170-
q_eff_fix_tmp = Array{T, 3}(undef, length(mpx), length(mpy), length(mpz))
171-
for iz in 1:size(ϵ, 3)
172-
pos_z::T = mpz[iz]
173-
for iy in 1:size(ϵ, 2)
174-
pos_y::T = mpy[iy]
175-
for ix in 1:size(ϵ, 1)
176-
pos_x = mpx[ix]
177-
pt::CartesianPoint{T} = CartesianPoint{T}(pos_x, pos_y, pos_z)
178-
ρ_tmp[ix, iy, iz]::T, ϵ[ix, iy, iz]::T, q_eff_fix_tmp[ix, iy, iz]::T = get_ρ_and_ϵ(pt, det, medium)
179-
end
178+
medium_ϵ_r::T = medium.ϵ_r
179+
ϵ = fill(medium_ϵ_r, length(mpx), length(mpy), length(mpz))
180+
ρ_tmp = zeros(T, length(mpx), length(mpy), length(mpz))
181+
q_eff_fix_tmp = zeros(T, length(mpx), length(mpy), length(mpz))
182+
fill_ρ_and_ϵ!(ϵ, ρ_tmp, q_eff_fix_tmp, Cartesian, mpz, mpy, mpx, use_nthreads, det.semiconductor)
183+
if !ismissing(det.passives)
184+
for passive in det.passives
185+
fill_ρ_and_ϵ!(ϵ, ρ_tmp, q_eff_fix_tmp, Cartesian, mpz, mpy, mpx, use_nthreads, passive)
180186
end
181187
end
182188
ϵ0_inv::T = inv(ϵ0)
@@ -272,8 +278,11 @@ function PotentialSimulationSetupRB(det::SolidStateDetector{T}, grid::CartesianG
272278

273279
potential::Array{T, 3} = ismissing(potential_array) ? zeros(T, size(grid)...) : potential_array
274280
point_types::Array{PointType, 3} = ones(PointType, size(grid)...)
275-
set_point_types_and_fixed_potentials!( point_types, potential, grid, det, weighting_potential_contact_id = weighting_potential_contact_id,
276-
not_only_paint_contacts = Val(not_only_paint_contacts), paint_contacts = Val(paint_contacts) )
281+
set_point_types_and_fixed_potentials!( point_types, potential, grid, det,
282+
weighting_potential_contact_id = weighting_potential_contact_id,
283+
use_nthreads = use_nthreads,
284+
not_only_paint_contacts = Val(not_only_paint_contacts),
285+
paint_contacts = Val(paint_contacts) )
277286
rbpotential::Array{T, 4} = RBExtBy2Array( potential, grid )
278287
rbpoint_types::Array{T, 4} = RBExtBy2Array( point_types, grid )
279288
potential = clear(potential)

0 commit comments

Comments
 (0)