Using pyplot@Plots on Julia@spack

 Using pyplot@Plots on Julia installed via spack is a bit complicated.

We first need to setup Python environment
(actually, you may want to use conda provided by Julia, so, you may do
Julia> ENV["PYTHON"]=""
)

as described in the following link;
https://github.com/JuliaPy/PyCall.jl

and then,

Julia> using Plots
Julia> pyplot()

offsetarray as a member of structure in Julia

 I do feel that accessing offsetarray, which is a member of a structure in Julia, is extremely slow.

In my case, I use a staggered grid in a FDTD code, and defined a structure as;

struct grid

    p :: OffsetArray{Float64,2}

    vx:: OffsetArray{Float64,2}

    vy:: OffsetArray{Float64,2}

end

then accessed like (output, input, resOfCalcF are all "grid".)

for j in -m+1:gridParam.ny+m

    for i in -m+1:gridParam.nx+m

         output.p[i,j] = input.p[i,j] + coef*resOfCalcF.p[i,j]

    end

end

It seems that Julia cannot determine the type of ".p" and this causes the type-instability.

To avoid this, I needed to cast type as 

op = parent(output.p)::Matrix{Float64}

ip = parent(input.p)::Matrix{Float64}

rp = parent(resOfCalcF.p)::Matrix{Float64}

and the array starts from 1 (the loop starts from 2, because the outer most layer is just a dummy)

for j in 1+1:gridParam.ny+m+m+1

    for i in 1+1:gridParam.nx+m+m+1

         op[i,j] = ip[i,j] + coef*rp[i,j] :: Float64

    end

end

I have not measured the speed-up, but probably x10 - x100 speedup is obtained.

FDM for the time dependent wave equation in Julia

I am writing a test code of our new time-dependent wave equation in Julia.

For now, it is equivalent to the normal wave equation and it looks fine.

I am not sure Julia is that fast and easy to write, but one of the merits is that the visualization functions can be called within Julia's framework.

 

Julia

この投稿が面白かったので、Juliaの練習を兼ねてちょっと遊んで見る。
https://twitter.com/solove_math/status/1035480802144833536

==ここからコード==
using Plots

function func2(n)
      return [1.0/(2.0*Float64(n) - 1.0),-1.0/(4.0*Float64(n) - 2.0) ,-1.0/(4.0*Float64(n)) ]
end

function func1(n)
      return [1.0/(2.0*Float64(n) - 1.0),-1.0/(2.0*Float64(n) )]
end

function main()

   n = 60
   y = 0.0
   y2= 0.0
   arr_y = zeros(Float64,n, 2)
   #arr_x = zeros(Int64,n)
   arr_x = 1:n
   for i in 1:n/2
      y = func1(i)
      arr_y[2*Int64(i)-1,1]= y[1]
      arr_y[2*Int64(i)  ,1]= y[2]
   end
   for i in 1:n/3
      y = func2(i)
      arr_y[3*Int64(i)-2,2]= y[1]
      arr_y[3*Int64(i)-1,2]= y[2]
      arr_y[3*Int64(i)  ,2]= y[3]
   end
   for i in 2:n
      arr_y[i,1]= arr_y[i,1] + arr_y[i-1,1]
      arr_y[i,2]= arr_y[i,2] + arr_y[i-1,2]
   end

   plt = plot(arr_x,arr_y)
   gui(plt)
   readline()
end


main()
==ここまでコード==