Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bugfix: Don't use symmetries for weighting potentials #216

Merged
merged 3 commits into from
Sep 14, 2021
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 13 additions & 11 deletions src/Simulation/Simulation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -239,7 +239,13 @@ function Grid(sim::Simulation{T, Cylindrical};

world_Δr, world_Δφ, world_Δz = width.(world.intervals)
world_r_mid = (world.intervals[1].right + world.intervals[1].left)/2

if for_weighting_potential && world_Δφ > 0
world_φ_int = SSDInterval{T, :closed, :open, :periodic, :periodic}(0, 2π)
world_Δφ = width(world_φ_int)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we always want it to be full ? I think this makes the full_2π keyword obsolete?

We could set full_2π = true here...

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I think it is obsolete. It's nowhere used. I will remove it.

I would keep it as we did it before. Always calculate the weigting potentials for 2π.
Only if it is completly phi-symmetric. Than, also only 2D.

I think we should refactor the symmetry handling, #215, in general. But for now we just need a quick fix.

else
world_φ_int = world.intervals[2]
end

max_distance_z = T(world_Δz / 4)
max_distance_φ = T(world_Δφ / 4)
max_distance_r = T(world_Δr / 4)
Expand Down Expand Up @@ -278,13 +284,13 @@ function Grid(sim::Simulation{T, Cylindrical};
important_z_points = initialize_axis_ticks(important_z_points; max_ratio = T(max_distance_ratio))
important_z_points = fill_up_ticks(important_z_points, max_distance_z)

append!(important_φ_points, endpoints(world.intervals[2])...)
append!(important_φ_points, endpoints(world_φ_int)...)
important_φ_points = unique!(sort!(important_φ_points))
if add_points_between_important_point
important_φ_points = sort!(vcat(important_φ_points, StatsBase.midpoints(important_φ_points)))
end
iL = searchsortedfirst(important_φ_points, world.intervals[2].left)
iR = searchsortedfirst(important_φ_points, world.intervals[2].right)
iL = searchsortedfirst(important_φ_points, world_φ_int.left)
iR = searchsortedfirst(important_φ_points, world_φ_int.right)
important_φ_points = unique(map(t -> isapprox(t, 0, atol = 1e-3) ? zero(T) : t, important_φ_points[iL:iR]))
important_φ_points = merge_close_ticks(important_φ_points, min_diff = T(1e-3))
important_φ_points = initialize_axis_ticks(important_φ_points; max_ratio = T(max_distance_ratio))
Expand All @@ -296,9 +302,9 @@ function Grid(sim::Simulation{T, Cylindrical};
ax_r = even_tick_axis(DiscreteAxis{T, BL, BR}(int_r, important_r_points))

# φ
L, R, BL, BR = get_boundary_types(world.intervals[2])
int_φ = Interval{L, R, T}(endpoints(world.intervals[2])...)
if full_2π || (for_weighting_potential && (world.intervals[2].left != world.intervals[2].right))
L, R, BL, BR = get_boundary_types(world_φ_int)
int_φ = Interval{L, R, T}(endpoints(world_φ_int)...)
if full_2π || (for_weighting_potential && (world_φ_int.left != world_φ_int.right))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

.. and then we can check only the statement if full_2π here?

L, R, BL, BR = :closed, :open, :periodic, :periodic
int_φ = Interval{L, R, T}(0, 2π)
end
Expand Down Expand Up @@ -1206,7 +1212,6 @@ There are several keyword arguments which can be used to tune the simulation.
* `max_distance_ratio::Real`: Maximum allowed ratio between the two distances in any dimension to the two neighbouring grid points.
If the ratio is too large, additional ticks are generated such that the new ratios are smaller than `max_distance_ratio`.
Default is `5`.
* `grid::Grid`: Initial grid used to start the simulation. Default is `Grid(sim)`.
* `depletion_handling::Bool`: Enables the handling of undepleted regions. Default is `false`.
* `use_nthreads::Int`: Number of threads to use in the computation. Default is `Base.Threads.nthreads()`.
The environment variable `JULIA_NUM_THREADS` must be set appropriately before the Julia session was
Expand Down Expand Up @@ -1238,7 +1243,6 @@ function simulate!( sim::Simulation{T, S};
min_tick_distance::Union{Missing, LengthQuantity, Tuple{LengthQuantity, AngleQuantity, LengthQuantity}} = missing,
max_tick_distance::Union{Missing, LengthQuantity, Tuple{LengthQuantity, AngleQuantity, LengthQuantity}} = missing,
max_distance_ratio::Real = 5,
grid::Grid{T, 3, S} = Grid(sim),
depletion_handling::Bool = false,
use_nthreads::Union{Int, Vector{Int}} = Base.Threads.nthreads(),
sor_consts::Union{Missing, <:Real, Tuple{<:Real,<:Real}} = missing,
Expand All @@ -1252,7 +1256,6 @@ function simulate!( sim::Simulation{T, S};
min_tick_distance = min_tick_distance,
max_tick_distance = max_tick_distance,
max_distance_ratio = max_distance_ratio,
grid = grid,
depletion_handling = depletion_handling,
use_nthreads = use_nthreads,
sor_consts = sor_consts,
Expand All @@ -1268,7 +1271,6 @@ function simulate!( sim::Simulation{T, S};
min_tick_distance = min_tick_distance,
max_tick_distance = max_tick_distance,
max_distance_ratio = max_distance_ratio,
grid = grid,
depletion_handling = depletion_handling,
use_nthreads = use_nthreads,
sor_consts = sor_consts,
Expand Down
4 changes: 2 additions & 2 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -118,9 +118,9 @@ end
@testset "Simulate example detector: Toroidal" begin
sim = Simulation{T}(SSD_examples[:CoaxialTorus])
SolidStateDetectors.apply_initial_state!(sim, ElectricPotential)
simulate!(sim, convergence_limit = 1e-6, refinement_limits = [0.2, 0.1, 0.05],
simulate!(sim, convergence_limit = 1e-5, refinement_limits = [0.2, 0.1, 0.05, 0.02, 0.01],
max_tick_distance = 0.5u"mm", verbose = false)
evt = Event([CartesianPoint{T}(0.01,0,0.003)])
evt = Event([CartesianPoint{T}(0.0075,0,0)])
simulate!(evt, sim, Δt = 1e-9, max_nsteps = 10000)
signalsum = T(0)
for i in 1:length(evt.waveforms)
Expand Down