11# 1D diffusion problem
22
33# Packages and inclusions
4- using ModelingToolkit,DiffEqOperators,LinearAlgebra,Test,OrdinaryDiffEq
4+ using ModelingToolkit,DiffEqOperators,LinearAlgebra,Test
55using ModelingToolkit: Differential
6- using DifferentialEquations
6+ using NonlinearSolve
77using DomainSets
88
99
@@ -17,7 +17,7 @@ using DomainSets
1717 ns = NonlinearSystem (eqs, [a, b, c, d], [])
1818 f = eval (generate_function (ns, [a, b, c, d])[2 ])
1919 prob = NonlinearProblem (ns, zeros (4 ), [])
20- sol = solve (prob)
20+ sol = NonlinearSolve . solve (prob, NewtonRaphson () )
2121 @test sol. u ≈ ones (4 )
2222end
2323
3939 pdesys = PDESystem ([eq],bcs,domains,[x],[u (x)])
4040 discretization = MOLFiniteDifference ([x=> dx], nothing , centered_order= 2 )
4141 prob = discretize (pdesys,discretization)
42- sol = solve (prob)
42+ sol = NonlinearSolve . solve (prob, NewtonRaphson () )
4343
4444 @test sol. u ≈ ones (4 )
4545end
6262 pdesys = PDESystem ([eq],bcs,domains,[x],[u (x)])
6363 discretization = MOLFiniteDifference ([x=> dx], nothing , centered_order= 2 )
6464 prob = discretize (pdesys,discretization)
65- sol = solve (prob)
65+ sol = NonlinearSolve . solve (prob, NewtonRaphson () )
6666
6767 @test sol. u ≈ 1.0 : 0.1 : 2.0
6868end
9595 discretization = MOLFiniteDifference ([x=> dx,y=> dy], nothing , centered_order= 2 )
9696
9797 prob = discretize (pdesys,discretization)
98- sol = solve (prob)
98+ sol = NonlinearSolve . solve (prob, NewtonRaphson () )
9999 xs,ys = [infimum (d. domain): dx: supremum (d. domain) for d in domains]
100100 u_sol = reshape (sol. u, (length (xs),length (ys)))
101101
102102 # test boundary
103- @test all (u_sol[:,1 ] .< eps ())
104- @test all (u_sol[1 ,:] .< eps ())
103+ @test all (abs .( u_sol[:,1 ]) .< eps (Float32 ))
104+ @test all (abs .( u_sol[1 ,:]) .< eps (Float32 ))
105105 @test u_sol[:,end ] ≈ 0 : dy: 1.0
106106 @test u_sol[end ,:] ≈ 0 : dx: 1.0
107107
108108 # test interior with finite differences
109109 interior = CartesianIndices ((axes (xs)[1 ], axes (ys)[1 ]))[2 : end - 1 ,2 : end - 1 ]
110110 fd = map (interior) do I
111- u_sol[(I - CartesianIndex (1 , 0 ))] + u_sol[(I + CartesianIndex (1 ,0 ))] - 2 * u_sol[I] < eps (Float32)
111+ abs ( u_sol[(I - CartesianIndex (1 , 0 ))] + u_sol[(I + CartesianIndex (1 ,0 ))] - 2 * u_sol[I]) < eps (Float32)
112112 end
113113 @test all (fd)
114114end
0 commit comments