|
| 1 | +using NonlinearSolve, LinearAlgebra, SparseArrays |
| 2 | + |
| 3 | +const N = 32 |
| 4 | +const xyd_brusselator = range(0,stop=1,length=N) |
| 5 | +brusselator_f(x, y) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * 5. |
| 6 | +limit(a, N) = a == N+1 ? 1 : a == 0 ? N : a |
| 7 | +function brusselator_2d_loop(du, u, p) |
| 8 | + A, B, alpha, dx = p |
| 9 | + alpha = alpha/dx^2 |
| 10 | + @inbounds for I in CartesianIndices((N, N)) |
| 11 | + i, j = Tuple(I) |
| 12 | + x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]] |
| 13 | + ip1, im1, jp1, jm1 = limit(i+1, N), limit(i-1, N), limit(j+1, N), limit(j-1, N) |
| 14 | + du[i,j,1] = alpha*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) + |
| 15 | + B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y) |
| 16 | + du[i,j,2] = alpha*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) + |
| 17 | + A*u[i,j,1] - u[i,j,1]^2*u[i,j,2] |
| 18 | + end |
| 19 | +end |
| 20 | +p = (3.4, 1., 10., step(xyd_brusselator)) |
| 21 | + |
| 22 | +function init_brusselator_2d(xyd) |
| 23 | + N = length(xyd) |
| 24 | + u = zeros(N, N, 2) |
| 25 | + for I in CartesianIndices((N, N)) |
| 26 | + x = xyd[I[1]] |
| 27 | + y = xyd[I[2]] |
| 28 | + u[I,1] = 22*(y*(1-y))^(3/2) |
| 29 | + u[I,2] = 27*(x*(1-x))^(3/2) |
| 30 | + end |
| 31 | + u |
| 32 | +end |
| 33 | +u0 = init_brusselator_2d(xyd_brusselator) |
| 34 | +prob_brusselator_2d = NonlinearProblem(brusselator_2d_loop,u0,p) |
| 35 | +sol = solve(prob_brusselator_2d, NewtonRaphson()) |
| 36 | + |
| 37 | +using Symbolics |
| 38 | +du0 = copy(u0) |
| 39 | +jac_sparsity = Symbolics.jacobian_sparsity((du,u)->brusselator_2d_loop(du,u,p),du0,u0) |
| 40 | + |
| 41 | +f = NonlinearFunction(brusselator_2d_loop;jac_prototype=float.(jac_sparsity)) |
| 42 | +prob_brusselator_2d = NonlinearProblem(f,u0,p) |
| 43 | +sol = solve(prob_brusselator_2d, NewtonRaphson()) |
0 commit comments