|
1 | 1 | using NonlinearSolve, LinearAlgebra, SparseArrays, Symbolics
|
2 | 2 |
|
3 | 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 |
| 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.0 |
| 6 | +limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a |
7 | 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] |
| 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), |
| 14 | + limit(j - 1, N) |
| 15 | + du[i, j, 1] = alpha * (u[im1, j, 1] + u[ip1, j, 1] + u[i, jp1, 1] + u[i, jm1, 1] - |
| 16 | + 4u[i, j, 1]) + |
| 17 | + B + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] + |
| 18 | + brusselator_f(x, y) |
| 19 | + du[i, j, 2] = alpha * (u[im1, j, 2] + u[ip1, j, 2] + u[i, jp1, 2] + u[i, jm1, 2] - |
| 20 | + 4u[i, j, 2]) + |
| 21 | + A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2] |
18 | 22 | end
|
19 | 23 | end
|
20 |
| -p = (3.4, 1., 10., step(xyd_brusselator)) |
| 24 | +p = (3.4, 1.0, 10.0, step(xyd_brusselator)) |
21 | 25 |
|
22 | 26 | 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 |
| 27 | + N = length(xyd) |
| 28 | + u = zeros(N, N, 2) |
| 29 | + for I in CartesianIndices((N, N)) |
| 30 | + x = xyd[I[1]] |
| 31 | + y = xyd[I[2]] |
| 32 | + u[I, 1] = 22 * (y * (1 - y))^(3 / 2) |
| 33 | + u[I, 2] = 27 * (x * (1 - x))^(3 / 2) |
| 34 | + end |
| 35 | + u |
32 | 36 | end
|
33 | 37 | u0 = init_brusselator_2d(xyd_brusselator)
|
34 |
| -prob_brusselator_2d = NonlinearProblem(brusselator_2d_loop,u0,p) |
| 38 | +prob_brusselator_2d = NonlinearProblem(brusselator_2d_loop, u0, p) |
35 | 39 | sol = solve(prob_brusselator_2d, NewtonRaphson())
|
36 | 40 |
|
37 | 41 | du0 = copy(u0)
|
38 |
| -jac_sparsity = Symbolics.jacobian_sparsity((du,u)->brusselator_2d_loop(du,u,p),du0,u0) |
| 42 | +jac_sparsity = Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p), du0, |
| 43 | + u0) |
39 | 44 |
|
40 |
| -ff = NonlinearFunction(brusselator_2d_loop;jac_prototype=float.(jac_sparsity)) |
41 |
| -prob_brusselator_2d = NonlinearProblem(ff,u0,p) |
| 45 | +ff = NonlinearFunction(brusselator_2d_loop; jac_prototype = float.(jac_sparsity)) |
| 46 | +prob_brusselator_2d = NonlinearProblem(ff, u0, p) |
42 | 47 | sol = solve(prob_brusselator_2d, NewtonRaphson())
|
43 | 48 | @test norm(sol.resid) < 1e-8
|
44 | 49 |
|
45 |
| -sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff=false)) |
| 50 | +sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = false)) |
46 | 51 | @test norm(sol.resid) < 1e-6
|
47 | 52 |
|
48 | 53 | cache = init(prob_brusselator_2d, NewtonRaphson())
|
|
0 commit comments