|
7 | 7 | @testitem "Manual SCC" setup=[CoreRootfindTesting] tags=[:core] begin
|
8 | 8 | using NonlinearSolveFirstOrder
|
9 | 9 | function f(du, u, p)
|
10 |
| - du[1] = cos(u[2]) - u[1] |
11 |
| - du[2] = sin(u[1] + u[2]) + u[2] |
12 |
| - du[3] = 2u[4] + u[3] + 1.0 |
13 |
| - du[4] = u[5]^2 + u[4] |
14 |
| - du[5] = u[3]^2 + u[5] |
15 |
| - du[6] = u[1] + u[2] + u[3] + u[4] + u[5] + 2.0u[6] + 2.5u[7] + 1.5u[8] |
16 |
| - du[7] = u[1] + u[2] + u[3] + 2.0u[4] + u[5] + 4.0u[6] - 1.5u[7] + 1.5u[8] |
17 |
| - du[8] = u[1] + 2.0u[2] + 3.0u[3] + 5.0u[4] + 6.0u[5] + u[6] - u[7] - u[8] |
| 10 | + du[1]=cos(u[2])-u[1] |
| 11 | + du[2]=sin(u[1]+u[2])+u[2] |
| 12 | + du[3]=2u[4]+u[3]+1.0 |
| 13 | + du[4]=u[5]^2+u[4] |
| 14 | + du[5]=u[3]^2+u[5] |
| 15 | + du[6]=u[1]+u[2]+u[3]+u[4]+u[5]+2.0u[6]+2.5u[7]+1.5u[8] |
| 16 | + du[7]=u[1]+u[2]+u[3]+2.0u[4]+u[5]+4.0u[6]-1.5u[7]+1.5u[8] |
| 17 | + du[8]=u[1]+2.0u[2]+3.0u[3]+5.0u[4]+6.0u[5]+u[6]-u[7]-u[8] |
18 | 18 | end
|
19 |
| - prob = NonlinearProblem(f, zeros(8)) |
20 |
| - sol = solve(prob, NewtonRaphson()) |
| 19 | + prob=NonlinearProblem(f, zeros(8)) |
| 20 | + sol=solve(prob, NewtonRaphson()) |
21 | 21 |
|
22 |
| - u0 = zeros(2) |
23 |
| - p = zeros(3) |
| 22 | + u0=zeros(2) |
| 23 | + p=zeros(3) |
24 | 24 |
|
25 | 25 | function f1(du, u, p)
|
26 |
| - du[1] = cos(u[2]) - u[1] |
27 |
| - du[2] = sin(u[1] + u[2]) + u[2] |
| 26 | + du[1]=cos(u[2])-u[1] |
| 27 | + du[2]=sin(u[1]+u[2])+u[2] |
28 | 28 | end
|
29 |
| - explicitfun1(p, sols) = nothing |
30 |
| - prob1 = NonlinearProblem( |
| 29 | + explicitfun1(p, sols)=nothing |
| 30 | + prob1=NonlinearProblem( |
31 | 31 | NonlinearFunction{true, SciMLBase.NoSpecialize}(f1), zeros(2), p)
|
32 |
| - sol1 = solve(prob1, NewtonRaphson()) |
| 32 | + sol1=solve(prob1, NewtonRaphson()) |
33 | 33 |
|
34 | 34 | function f2(du, u, p)
|
35 |
| - du[1] = 2u[2] + u[1] + 1.0 |
36 |
| - du[2] = u[3]^2 + u[2] |
37 |
| - du[3] = u[1]^2 + u[3] |
| 35 | + du[1]=2u[2]+u[1]+1.0 |
| 36 | + du[2]=u[3]^2+u[2] |
| 37 | + du[3]=u[1]^2+u[3] |
38 | 38 | end
|
39 |
| - explicitfun2(p, sols) = nothing |
40 |
| - prob2 = NonlinearProblem( |
| 39 | + explicitfun2(p, sols)=nothing |
| 40 | + prob2=NonlinearProblem( |
41 | 41 | NonlinearFunction{true, SciMLBase.NoSpecialize}(f2), zeros(3), p)
|
42 |
| - sol2 = solve(prob2, NewtonRaphson()) |
| 42 | + sol2=solve(prob2, NewtonRaphson()) |
43 | 43 |
|
44 | 44 | # Convert f3 to a LinearProblem since it's linear in u
|
45 | 45 | # du = Au + b where A is the coefficient matrix and b is from parameters
|
46 |
| - A3 = [2.0 2.5 1.5; 4.0 -1.5 1.5; 1.0 -1.0 -1.0] |
47 |
| - b3 = p # b will be updated by explicitfun3 |
48 |
| - prob3 = LinearProblem(A3, b3, zeros(3)) |
| 46 | + A3=[2.0 2.5 1.5; 4.0 -1.5 1.5; 1.0 -1.0 -1.0] |
| 47 | + b3=p # b will be updated by explicitfun3 |
| 48 | + prob3=LinearProblem(A3, b3, zeros(3)) |
49 | 49 | function explicitfun3(p, sols)
|
50 |
| - p[1] = -(sols[1][1] + sols[1][2] + sols[2][1] + sols[2][2] + sols[2][3]) |
51 |
| - p[2] = -(sols[1][1] + sols[1][2] + sols[2][1] + 2.0sols[2][2] + sols[2][3]) |
52 |
| - p[3] = -(sols[1][1] + 2.0sols[1][2] + 3.0sols[2][1] + 5.0sols[2][2] + |
| 50 | + p[1]=-(sols[1][1]+sols[1][2]+sols[2][1]+sols[2][2]+sols[2][3]) |
| 51 | + p[2]=-(sols[1][1]+sols[1][2]+sols[2][1]+2.0sols[2][2]+sols[2][3]) |
| 52 | + p[3]=-(sols[1][1]+2.0sols[1][2]+3.0sols[2][1]+5.0sols[2][2]+ |
53 | 53 | 6.0sols[2][3])
|
54 | 54 | end
|
55 | 55 | explicitfun3(p, [sol1, sol2])
|
56 |
| - sol3 = solve(prob3) # LinearProblem uses default linear solver |
57 |
| - manualscc = reduce(vcat,(sol1, sol2, sol3)) |
| 56 | + sol3=solve(prob3) # LinearProblem uses default linear solver |
| 57 | + manualscc=reduce(vcat, (sol1, sol2, sol3)) |
58 | 58 |
|
59 |
| - sccprob = SciMLBase.SCCNonlinearProblem((prob1, prob2, prob3), |
| 59 | + sccprob=SciMLBase.SCCNonlinearProblem((prob1, prob2, prob3), |
60 | 60 | SciMLBase.Void{Any}.([explicitfun1, explicitfun2, explicitfun3]))
|
61 |
| - |
| 61 | + |
62 | 62 | # Test with SCCAlg that handles both nonlinear and linear problems
|
63 | 63 | using SCCNonlinearSolve
|
64 |
| - scc_alg = SCCNonlinearSolve.SCCAlg(nlalg = NewtonRaphson(), linalg = nothing) |
65 |
| - scc_sol = solve(sccprob, scc_alg) |
| 64 | + scc_alg=SCCNonlinearSolve.SCCAlg(nlalg = NewtonRaphson(), linalg = nothing) |
| 65 | + scc_sol=solve(sccprob, scc_alg) |
66 | 66 | @test sol ≈ manualscc ≈ scc_sol
|
67 | 67 |
|
68 | 68 | import NonlinearSolve # Required for Default
|
69 | 69 |
|
70 | 70 | # Test default interface
|
71 |
| - scc_sol_default = solve(sccprob) |
| 71 | + scc_sol_default=solve(sccprob) |
72 | 72 | @test sol ≈ manualscc ≈ scc_sol_default
|
73 | 73 | end
|
0 commit comments