julia - 交互张量的积分
问题描述
我需要在每个时间步计算以下积分数千次:
到目前为止,我已经在 Julia 中实现了:
using StaticArrays
function interactiontensor(C, a1, a2, a3, ϕ, θ)
n1,n2 = 100,50
T = fill(0.0,3,3,3,3)
Av = zeros(4,4)
invAv = similar(Av)
xi = Vector{Float64}(undef, 3)
@inbounds for p ∈ 1:n1
sinθp = sind(θ[p])
cosθp = cosd(θ[p])
for q ∈ 1:n2
sinϕq = sind(ϕ[q])
cosϕq = cosd(ϕ[q])
# -- Director cosines
xi[1] = sinθp*cosϕq/a1
xi[2] = sinθp*sinϕq/a2
xi[3] = cosθp/a3
Christoffel!(Av,C,xi)
fillAv!(Av, xi)
invAv = inv(SMatrix{4,4}(Av))
tensorT!(T,invAv,xi,sinθp)
surface += sinθp
end
end
return T ./= surface
end
@inline function Christoffel!(Av,C,xi)
@inbounds for t ∈ 1:3, r ∈ 1:3
aux = zero(eltype(C))
for u ∈ 1:3, s ∈ 1:3
aux += C[r, s, t, u] * xi[s] * xi[u]
end
Av[r, t] = aux
end
end
@inline function tensorT!(T,invAv,xi,sinθp)
@inbounds for k ∈ 1:3, i ∈ 1:3
aux = invAv[i, k]
for l ∈ 1:3, j ∈ 1:3
T[i, j, k, l] += aux * xi[j] * xi[l] * sinθp
end
end
end
@inline function fillAv!(Av, xi)
@inbounds for i ∈ 1:3
xi0 = xi[i]
Av[i, 4] = xi0
Av[4, i] = xi0
end
end
和
n1,n2 = 100,100
step = π/n1
dθ,dϕ = π/n1, 2π/n2
θ = rad2deg.(range(dθ, stop = pi, length = n1))
ϕ = rad2deg.(range(dϕ, stop = 2pi, length = n2))
C = @SArray rand(3,3,3,3)
@btime interactiontensor($C, $10.0, $5.0, $1.0, $ϕ, $θ);
# 544.795 μs (4 allocations: 1.08 KiB)
考虑到我理想地需要计算这个积分的次数,我的实现是否有任何优化或替代方法,以显着降低计算成本?
解决方案
以下是一些建议:
sinθp, cosθp = sincosd(θ[p])
,即一步计算正弦和余弦。- 初始化
xi = @SVector zeros(3)
为静态向量,然后使用Setfield.jl在每次迭代中分配值,即@set x[1] = sinθp*cosϕq/a1
. - 加载包LoopVectorization.jl并使用
@avx
宏(非常粗略地说类似于@simd
)来加速和中Christoffel!
的循环。tensorT!
fillAv!
在我的机器上,我发现这些更改将计算时间减少了 5 倍以上(相对于 OP 中的原始函数)。最大的一块是由于@avx
,上面的第二点大约是~30%。
julia> @btime interactiontensor_original($C, $10.0, $5.0, $1.0, $ϕ, $θ);
661.655 μs (5 allocations: 1.28 KiB)
julia> @btime interactiontensor_optimized($C, $10.0, $5.0, $1.0, $ϕ, $θ);
125.352 μs (4 allocations: 1.17 KiB)
这里是完整的修改代码(请注意,我注释掉了surface
OP 中未指定的行):
using StaticArrays, Setfield, LoopVectorization
function interactiontensor_optimized(C, a1, a2, a3, ϕ, θ)
n1,n2 = 100,50
T = fill(0.0,3,3,3,3)
Av = zeros(4,4)
invAv = similar(Av)
xi = @SVector zeros(3)
@inbounds for p ∈ 1:n1
sinθp, cosθp = sincosd(θ[p])
for q ∈ 1:n2
sinϕq, cosϕq = sincosd(ϕ[q])
# -- Director cosines
@set xi[1] = sinθp*cosϕq/a1
@set xi[2] = sinθp*sinϕq/a2
@set xi[3] = cosθp/a3
Christoffel!(Av,C,xi)
fillAv!(Av, xi)
invAv = inv(SMatrix{4,4}(Av))
tensorT!(T,invAv,xi,sinθp)
# surface += sinθp
end
end
return T #./= surface
end
@inline function Christoffel!(Av,C,xi)
@avx for t ∈ 1:3, r ∈ 1:3
aux = zero(eltype(C))
for u ∈ 1:3, s ∈ 1:3
aux += C[r, s, t, u] * xi[s] * xi[u]
end
Av[r, t] = aux
end
end
@inline function tensorT!(T,invAv,xi,sinθp)
@avx for k ∈ 1:3, i ∈ 1:3
aux = invAv[i, k]
for l ∈ 1:3, j ∈ 1:3
T[i, j, k, l] += aux * xi[j] * xi[l] * sinθp
end
end
end
@inline function fillAv!(Av, xi)
@avx for i ∈ 1:3
xi0 = xi[i]
Av[i, 4] = xi0
Av[4, i] = xi0
end
end
推荐阅读
- python - 无法理解此代码的输出:python list to 1-d numpy array
- javascript - Javascript均衡数组中的元素
- node.js - 如何根据猫鼬上的子填充属性查找对象
- c++ - 如何更改 sizeof(std::string)
- windows - Flutter 在 Windows 机器上使用 LF 作为行尾
- treesitter - 是否有用于解析任意长度列表的标准treesitter 构造?
- azure - 如何在 Azure 中为存储容器创建容量警报
- sql - 从基于表的查询中排除数据对?
- android - 在 Flutter 中使用 flutter_filereader 将中文标题更改为英文
- javascript - Flex Box,对齐最后 2 个项目,位于 flex 框的底部