Skip to content

Commit

Permalink
Merge pull request #46 from lmh91/analytic_solutions
Browse files Browse the repository at this point in the history 
Performance fix + Extend tests (comparison to analytic solutions)
  • Loading branch information
@lmh91
lmh91 authored May 16, 2020
2 parents 4624176 + 41dcfdc commit 05633f2
Show file tree
Hide file tree
Showing 13 changed files with 364 additions and 167 deletions.
132 changes: 132 additions & 0 deletions examples/example_detector_config_files/infinite_coaxial_capacitor_cartesian_coords.json
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
{
"name": "Infinite Coaxial Capacitor",
"units": {
"length": "mm",
"angle": "deg",
"potential": "V",
"temperature": "K"
},
"grid": {
"coordinates": "cartesian",
"axes": {
"x": {
"from": -35.0,
"to": 35.0,
"boundaries": "reflecting"
},
"y": {
"from": -35.0,
"to": 35.0,
"boundaries": "reflecting"
},
"z": {
"from": 0,
"to": 40.0,
"boundaries": "reflecting"
}
}
},
"medium": "HPGe",
"objects": [
{
"type": "semiconductor",
"material": "HPGe",
"bulk_type": "p",
"temperature": 77.0,
"charge_density_model": {
"name": "linear",
"r": {
"init": 0.0,
"gradient": 0.0
},
"phi": {
"init": 0.0,
"gradient": 0.0
},
"z": {
"init": 0.0,
"gradient": 0.0
}
},
"geometry": {
"type": "difference",
"parts": [
{
"name": "Initial Cylinder",
"type": "tube",
"r": {
"from": 0.0,
"to": 35.0
},
"phi": {
"from": 0.0,
"to": 360.0
},
"z": {
"from": -1.0,
"to": 41.0
}
},
{
"name": "borehole",
"type": "tube",
"r": {
"from": 0.0,
"to": 5.0
},
"phi": {
"from": 0.0,
"to": 360.0
},
"z": {
"from": -1.0,
"to": 41.0
}
}
]
}
},
{
"type": "contact",
"material": "HPGe",
"channel": 1,
"potential": 0.0,
"geometry": {
"type": "tube",
"r": {
"from": 5.0,
"to": 5.0
},
"phi": {
"from": 0.0,
"to": 360.0
},
"z": {
"from": -1.0,
"to": 41.0
}
}
},
{
"type": "contact",
"material": "HPGe",
"channel": 2,
"potential": 10.0,
"geometry": {
"type": "tube",
"r": {
"from": 35.0,
"to": 35.0
},
"phi": {
"from": 0.0,
"to": 360.0
},
"z": {
"from": -1.0,
"to": 41.0
}
}
}
]
}
52 changes: 17 additions & 35 deletions src/ElectricField/ElectricField.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,11 +56,11 @@ end

function get_electric_field_from_potential(ep::ElectricPotential{T, 3, :cylindrical}, pointtypes::PointTypes{T}, fieldvector_coordinates=:xyz)::ElectricField{T, 3, :cylindrical} where {T <: SSDFloat}
p = ep.data
axr::Vector{T} = collect(ep.grid.r)
axφ::Vector{T} = collect(ep.grid.φ)
axz::Vector{T} = collect(ep.grid.z)
axr::Vector{T} = collect(ep.grid.axes[1])
axφ::Vector{T} = collect(ep.grid.axes[2])
axz::Vector{T} = collect(ep.grid.axes[3])

cyclic::T = ep.grid.φ.interval.right
cyclic::T = ep.grid.axes[2].interval.right
ef = Array{SVector{3, T}}(undef, size(p)...)
for iz in 1:size(ef, 3)
for iφ in 1:size(ef, 2)
Expand Down Expand Up @@ -185,55 +185,37 @@ function convert_field_vectors_to_xyz(field::Array{SArray{Tuple{3},T,1,3},3}, φ
for iz in axes(field, 3)
for ir in axes(field, 1)
field_xyz[ir,iφ,iz] = Rα * field[ir,iφ,iz]
# field_xyz[ir,iφ,iz] = get_xyz_vector_from_field_vector(field[ir,iφ,iz], φ)
end
end
end
return field_xyz
end
#
# function get_xyz_vector_from_field_vector(vector, α)
# Rα = @SArray([cos(α) -1*sin(α) 0;sin(α) cos(α) 0;0 0 1])
# result = Rα*vector
# result
# end


function interpolated_scalarfield(ep::ScalarPotential{T, 3, :cylindrical}) where {T}
knots = ep.grid.axes[1].ticks, cat(ep.grid.axes[2].ticks,T(2π),dims=1), ep.grid.axes[3].ticks#(grid.r, grid.φ, grid.z)
@inbounds knots = ep.grid.axes[1].ticks, cat(ep.grid.axes[2].ticks,T(2π),dims=1), ep.grid.axes[3].ticks
ext_data = cat(ep.data, ep.data[:,1:1,:], dims=2)
i = interpolate(knots, ext_data, Gridded(Linear()))
vector_field_itp = extrapolate(i, (Interpolations.Line(), Periodic(), Interpolations.Line()))
return vector_field_itp
end
function interpolated_scalarfield(ep::ScalarPotential{T, 3, :cartesian}) where {T}
knots = ep.grid.axes#(grid.x, grid.y, grid.z)
@inbounds knots = ep.grid.axes[1].ticks, ep.grid.axes[2].ticks, ep.grid.axes[3].ticks
i = interpolate(knots, ep.data, Gridded(Linear()))
vector_field_itp = extrapolate(i, (Interpolations.Line(), Interpolations.Line(), Interpolations.Line()))
return vector_field_itp
end
#
# function interpolated_vectorfield(vectorfield::AbstractArray{<:SVector{3, T},3}, ep::ElectricPotential{T}) where {T}
# knots = ep.grid.axes#(grid.r, grid.φ, grid.z)
# i = interpolate(knots, vectorfield, Gridded(Linear()))
# vectorfield_itp = extrapolate(i, Periodic())
# return vectorfield_itp
# end
#
# function setup_interpolated_vectorfield(vectorfield, grid::CylindricalGrid{T}) where {T}
# knots = grid.axes #(grid.r, grid.φ, grid.z)
# i = interpolate(knots, vectorfield, Gridded(Linear()))
# vector_field_itp = extrapolate(i, Periodic())
# return vector_field_itp
# end


function get_interpolated_drift_field(velocityfield, grid::CylindricalGrid{T}) where {T}
extended_velocityfield = cat(velocityfield, velocityfield[:,1:1,:], dims=2)
knots = grid.axes[1].ticks, cat(grid.axes[2].ticks,T(2π),dims=1), grid.axes[3].ticks
@inbounds knots = grid.axes[1].ticks, cat(grid.axes[2].ticks,T(2π),dims=1), grid.axes[3].ticks
i = interpolate(knots, extended_velocityfield, Gridded(Linear()))
velocity_field_itp = extrapolate(i, (Interpolations.Line(), Periodic(), Interpolations.Line()))
return velocity_field_itp
end
function get_interpolated_drift_field(velocityfield, grid::CartesianGrid{T}) where {T}
knots = grid.axes[:1].ticks, grid.axes[:2].ticks, grid.axes[:3].ticks
@inbounds knots = grid.axes[1].ticks, grid.axes[2].ticks, grid.axes[3].ticks
i = interpolate(knots, velocityfield, Gridded(Linear()))
velocity_field_itp = extrapolate(i, (Interpolations.Line(), Interpolations.Line(), Interpolations.Line()))
return velocity_field_itp
Expand All @@ -242,12 +224,12 @@ end
include("plot_recipes.jl")

function get_electric_field_from_potential(ep::ElectricPotential{T, 3, :cartesian}, pointtypes::PointTypes{T})::ElectricField{T, 3, :cartesian} where {T <: SSDFloat}
axx::Vector{T} = collect(ep.grid[:x])
axy::Vector{T} = collect(ep.grid[:y])
axz::Vector{T} = collect(ep.grid.z)
axx_ext::Vector{T} = get_extended_ticks(ep.grid[:x])
axy_ext::Vector{T} = get_extended_ticks(ep.grid[:y])
axz_ext::Vector{T} = get_extended_ticks(ep.grid.z)
axx::Vector{T} = collect(ep.grid.axes[1])
axy::Vector{T} = collect(ep.grid.axes[2])
axz::Vector{T} = collect(ep.grid.axes[3])
axx_ext::Vector{T} = get_extended_ticks(ep.grid.axes[1])
axy_ext::Vector{T} = get_extended_ticks(ep.grid.axes[2])
axz_ext::Vector{T} = get_extended_ticks(ep.grid.axes[3])

ef::Array{SVector{3, T}} = Array{SVector{3, T}}(undef, size(ep.data))

Expand Down
28 changes: 14 additions & 14 deletions src/ElectricField/plot_recipes.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,36 +32,36 @@ end
subplot := 1
title := "r_component"
ylabel --> "z ["*units[SI_factor]*"]"
grid.r ./ SI_factor, grid.z ./ SI_factor, vectorfield_r[:,i_fixed,:]'
grid.axes[1] ./ SI_factor, grid.axes[3] ./ SI_factor, vectorfield_r[:,i_fixed,:]'
end
@series begin
subplot := 2
colorbar_title --> "Field Strength [V / "*units[SI_factor]*"]"
title := "φ_component"

grid.r ./SI_factor, grid.z ./SI_factor, vectorfield_φ[:,i_fixed,:]'
grid.axes[1] ./SI_factor, grid.axes[3] ./SI_factor, vectorfield_φ[:,i_fixed,:]'
end
@series begin
subplot := 3
title := "z_component"
ylabel --> "z ["*units[SI_factor]*"]"
xlabel --> "r ["*units[SI_factor]*"]"
grid.r ./SI_factor, grid.z ./SI_factor, vectorfield_z[:,i_fixed,:]'
grid.axes[1] ./SI_factor, grid.axes[3] ./SI_factor, vectorfield_z[:,i_fixed,:]'
end
@series begin
subplot := 4
colorbar_title --> "Field Strength [V / "*units[SI_factor]*"]"
xlabel --> "r ["*units[SI_factor]*"]"
title:= "magnitude"
grid.r ./SI_factor, grid.z ./SI_factor, vectorfield_magn[:,i_fixed,:]'
grid.axes[1] ./SI_factor, grid.axes[3] ./SI_factor, vectorfield_magn[:,i_fixed,:]'
end
end
elseif view == :ef
if plane == :rφ
vectorfield_xyz = Array{Vector{T}}(undef,size(vectorfield,1),size(vectorfield,2),size(vectorfield,3));
for (iz,z) in enumerate(grid.z)
for (iφ,φ) in enumerate(grid.φ)
for (ir,r) in enumerate(grid.r)
for (iz,z) in enumerate(grid.axes[3])
for (iφ,φ) in enumerate(grid.axes[2])
for (ir,r) in enumerate(grid.axes[1])
vectorfield_xyz[ir,iφ,iz]=get_xyz_vector_from_rφz_field_vector_at_rφz(vectorfield,r,φ,z,ir,iφ,iz)
end
end
Expand All @@ -78,15 +78,15 @@ end
label := ""
ylabel := "y "
xlabel := "x "
title := "z = $(round(grid.z[i_fixed]/SI_factor,digits=2)) / mm"
xlims := (-1.2/SI_factor*maximum(grid.r),1.2/SI_factor*maximum(grid.r))
ylims := (-1.2/SI_factor*maximum(grid.r),1.2/SI_factor*maximum(grid.r))
for (ir,r) in enumerate(grid.r[1:spacing:end])
for (iφ,φ) in enumerate(grid.φ)
title := "z = $(round(grid.axes[3][i_fixed]/SI_factor,digits=2)) / mm"
xlims := (-1.2/SI_factor*maximum(grid.axes[1]),1.2/SI_factor*maximum(grid.axes[1]))
ylims := (-1.2/SI_factor*maximum(grid.axes[1]),1.2/SI_factor*maximum(grid.axes[1]))
for (ir,r) in enumerate(grid.axes[1][1:spacing:end])
for (iφ,φ) in enumerate(grid..axes[2])
x= r*cos(φ)
y= r*sin(φ)
ir_actual=findfirst(x->x==r,grid.r)
iφ_actual=findfirst(x->x==φ,grid.φ)
ir_actual=findfirst(x->x==r,grid.axes[1])
iφ_actual=findfirst(x->x==φ,grid.axes[2])
xy_magn = vectorfield_xy_magn[ir_actual,iφ_actual]
vector=vectorfield_xyz[ir_actual,iφ_actual,i_fixed]/xy_magn
vector*=((vectorfield_xy_magn[ir_actual,iφ_actual]-0.8*minimum(vectorfield_xy_magn))/diff_magn)
Expand Down
15 changes: 7 additions & 8 deletions src/Grids/Grids.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,8 @@ const SphericalGrid{T} = Grid{T, 3, :Spherical}
@inline getindex(g::Grid{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = broadcast(getindex, g.axes, I)
@inline getindex(g::Grid{T, N, S}, i::Int) where {T, N, S} = getproperty(g, :axes)[i]
@inline getindex(g::Grid{T, N, S}, s::Symbol) where {T, N, S} = getindex(g, Val{s}())
@inline getproperty(g::Grid{T, N, S}, s::Symbol) where {T, N, S} = getproperty(g, Val{s}())

@inline getproperty(g::Grid{T, N, S}, s::Symbol) where {T, N, S} = getproperty(g, Val{s}())
@inline getproperty(g::Grid{T}, ::Val{:axes}) where {T} = getfield(g, :axes)

@inline getproperty(g::CylindricalGrid{T}, ::Val{:axes}) where {T} = getfield(g, :axes)
Expand All @@ -35,7 +35,6 @@ const SphericalGrid{T} = Grid{T, 3, :Spherical}
@inline getproperty(g::CartesianGrid{T}, ::Val{:y}) where {T} = @inbounds g.axes[2]
@inline getproperty(g::CartesianGrid{T}, ::Val{:z}) where {T} = @inbounds g.axes[3]


@inline getindex(g::CylindricalGrid{T}, ::Val{:r}) where {T} = @inbounds g.axes[1]
@inline getindex(g::CylindricalGrid{T}, ::Val{:φ}) where {T} = @inbounds g.axes[2]
@inline getindex(g::CylindricalGrid{T}, ::Val{:z}) where {T} = @inbounds g.axes[3]
Expand Down Expand Up @@ -148,9 +147,9 @@ end
Base.convert(T::Type{Grid}, x::NamedTuple) = T(x)

function NamedTuple(grid::Grid{T, 3, :cylindrical}) where {T}
axr::DiscreteAxis{T} = grid.r
axφ::DiscreteAxis{T} = grid.φ
axz::DiscreteAxis{T} = grid.z
axr::DiscreteAxis{T} = grid.axes[1]
axφ::DiscreteAxis{T} = grid.axes[2]
axz::DiscreteAxis{T} = grid.axes[3]
return (
coordtype = "cylindrical",
ndims = 3,
Expand All @@ -162,9 +161,9 @@ function NamedTuple(grid::Grid{T, 3, :cylindrical}) where {T}
)
end
function NamedTuple(grid::Grid{T, 3, :cartesian}) where {T}
axx::DiscreteAxis{T} = grid[:x]
axy::DiscreteAxis{T} = grid[:y]
axz::DiscreteAxis{T} = grid.z
axx::DiscreteAxis{T} = grid.axes[1]
axy::DiscreteAxis{T} = grid.axes[2]
axz::DiscreteAxis{T} = grid.axes[3]
return (
coordtype = "cartesian",
ndims = 3,
Expand Down
Loading

0 comments on commit 05633f2

Please sign in to comment.