Skip to content

Vectorized Bruss #34

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

Open
YingboMa opened this issue Dec 2, 2018 · 1 comment
Open

Vectorized Bruss #34

YingboMa opened this issue Dec 2, 2018 · 1 comment
Assignees

Comments

@YingboMa
Copy link
Member

YingboMa commented Dec 2, 2018

 kernel_u! = let N=N, xyd=xyd_brusselator, dx=step(xyd_brusselator)
    @inline function (du, u, A, B, α, II, I, t)
      i, j = Tuple(I)
      x = xyd[I[1]]
      y = xyd[I[2]]
      ip1 = limit(i+1, N); im1 = limit(i-1, N)
      jp1 = limit(j+1, N); jm1 = limit(j-1, N)
      du[II[i,j,1]] = α[II[i,j,1]]*(u[II[im1,j,1]] + u[II[ip1,j,1]] + u[II[i,jp1,1]] + u[II[i,jm1,1]] - 4u[II[i,j,1]])/dx^2 +
      B[II[i,j,1]] + u[II[i,j,1]]^2*u[II[i,j,2]] - (A[II[i,j,1]] + 1)*u[II[i,j,1]] + brusselator_f(x, y, t)
    end
  end
  kernel_v! = let N=N, xyd=xyd_brusselator, dx=step(xyd_brusselator)
    @inline function (du, u, A, B, α, II, I, t)
      i, j = Tuple(I)
      ip1 = limit(i+1, N)
      im1 = limit(i-1, N)
      jp1 = limit(j+1, N)
      jm1 = limit(j-1, N)
      du[II[i,j,2]] = α[II[i,j,2]]*(u[II[im1,j,2]] + u[II[ip1,j,2]] + u[II[i,jp1,2]] + u[II[i,jm1,2]] - 4u[II[i,j,2]])/dx^2 +
      A[II[i,j,1]]*u[II[i,j,1]] - u[II[i,j,1]]^2*u[II[i,j,2]]
    end
  end
  brusselator_2d = let N=N, xyd=xyd_brusselator, dx=step(xyd_brusselator)
    function (du, u, p, t)
      @inbounds begin
        ii1 = N^2
        ii2 = ii1+N^2
        ii3 = ii2+2(N^2)
        A = @view p[1:ii1]
        B = @view p[ii1+1:ii2]
        α = @view p[ii2+1:ii3]
        II = LinearIndices((N, N, 2))
        kernel_u!.(Ref(du), Ref(u), Ref(A), Ref(B), Ref(α), Ref(II), CartesianIndices((N, N)), t)
        kernel_v!.(Ref(du), Ref(u), Ref(A), Ref(B), Ref(α), Ref(II), CartesianIndices((N, N)), t)
        return nothing
      end
    end
  end
@YingboMa YingboMa self-assigned this Dec 2, 2018
@ChrisRackauckas
Copy link
Member

@YingboMa can you PR this one today as a separate version of Bruss in here? This will be good for GPU tests I think.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants