# 加权Jacobi的迭代优化

function solvedByItSolvers!(::weightedJaboci{dim}, uh::Vector{Float64}, fh::Vector{Float64}, Ninter::Int, weight::Float64 = 0.5) where {dim}
uh = reshape(uh, Ninter+1, Ninter+1)
fh = reshape(fh, Ninter+1, Ninter+1)

# update Jacobi iteration
uhJ = uh
for ii = 2:Ninter
@simd for jj = 2:Ninter
@inbounds uhJ[ii,jj] = ( fh[ii,jj] + uh[ii-1,jj] + uh[ii+1,jj] + uh[ii,jj-1] + uh[ii,jj+1] )/4
end
end

uh = weight.*uh .+ (1-weight).*uhJ
return uh[:]
end



emmmm，简洁不代表快。因为Julia不能告诉LLVM数组的alias信息，目前。

using LinearAlgebra
using PyPlot
len = 300
x = y = range(-2pi,stop=2pi,length=len)
uh = [sin(xx^2)+cos(yy)^2 for xx in x, yy in y]
fh = [sin(xx)+cos(yy) for xx in x, yy in y]
function solve_poisson!(uh, fh, iter=1000)
@inbounds begin
n = LinearAlgebra.checksquare(uh)
II = LinearIndices((n, n))
for _ in 1:iter, ii in 2:n-1, jj in 2:n-1
uh[II[ii,jj]] = ( fh[II[ii,jj]] + uh[II[ii-1,jj]] + uh[II[ii+1,jj]] + uh[II[ii,jj-1]] + uh[II[ii,jj+1]] )/4
end
return uh
end
end
plt[:subplot](2, 1, 1)
contour(x, y, uh)
plt[:title]("Initial Condition")
solve_poisson!(uh, fh)
plt[:subplot](2, 1, 2)
contour(x, y, uh)
plt[:title]("After 1000 Gauss-Seidel Iterations")
plt[:savefig]("Poisson.png", dpi=600)


2 个赞

uhJ[ii,jj] = ( fh[ii,jj] + uh[ii-1,jj] + uh[ii+1,jj] + uh[ii,jj-1] + uh[ii,jj+1] )/4


for ii in 2:n-1, jj in 2:n-1
#---
end


for ii in 2:n-1
@simd for jj in 2:n-1
#---
end
end


ii = 2:n-1
jj = 2:n-1
uhJ[ii,jj] = ( fh[ii,jj] + uh[ii-1,jj] + uh[ii+1,jj] + uh[ii,jj-1] + uh[ii,jj+1] )/4


julia 不用向量化，但是 for 循环的写法上还挺多的，一般是怎么确定哪一种比较好呢？根据自己的经验？

julia> uh = zeros(5);

julia> uhJ = uh;

julia> pointer(uhJ)
Ptr{Float64} @0x00007fdd153fe7d0

julia> pointer(uh)
Ptr{Float64} @0x00007fdd153fe7d0

julia> uh[1] = 1000
1000

julia> uhJ
5-element Array{Float64,1}:
1000.0
0.0
0.0
0.0
0.0


for ii in 2:n-1, jj in 2:n-1
#---
end


for ii in 2:n-1
for jj in 2:n-1
# ---
end
end


2 个赞

using LinearAlgebra
using PyPlot
len = 300
x = y = range(-2pi,stop=2pi,length=len)
uh = [sin(xx^2)+cos(yy)^2 for xx in x, yy in y]
fh = [sin(xx)+cos(yy) for xx in x, yy in y]
function solve_poisson(uh, fh, iter=1000, weight=0.5)
@inbounds begin
n = LinearAlgebra.checksquare(uh)
uhJ = zero(uh)
II = LinearIndices((n, n))
for _ in 1:iter, ii in 2:n-1, jj in 2:n-1
uhJ[II[ii,jj]] = ( fh[II[ii,jj]] + uh[II[ii-1,jj]] + uh[II[ii+1,jj]] + uh[II[ii,jj-1]] + uh[II[ii,jj+1]] )/4
end
for ii in eachindex(uhJ)
uhJ[ii] = weight*uh[ii] + (1-weight)*uhJ[ii]
end
return uhJ
end
end
plt[:subplot](2, 1, 1)
contour(x, y, uh)
plt[:title]("Initial Condition")
uhJ = solve_poisson(uh, fh)
plt[:subplot](2, 1, 2)
contour(x, y, uhJ)
plt[:title]("After 1000 weighted Jacobi Iterations")
plt[:savefig]("Poisson.png", dpi=400)


1 个赞

DG 求解 Poisson？ 我看了 JuAFEM 的包，打算后面写个 DG 的