From 56bc1d9b565a80aad3e59833db08c99f0b48c66d Mon Sep 17 00:00:00 2001 From: romain veltz Date: Sun, 15 Sep 2024 20:49:55 +0200 Subject: [PATCH] !!!!! newton -> solve(, ::Newton,) --- examples/SH2d-fronts-cuda.jl | 4 ++-- examples/SH2d-fronts.jl | 6 +++--- examples/SH3d.jl | 2 +- examples/SHpde_snaking.jl | 2 +- examples/cGL2d-shooting.jl | 2 +- examples/cGL2d.jl | 2 +- examples/carrier.jl | 4 ++-- examples/chan-af.jl | 2 +- examples/chan.jl | 12 ++++++------ examples/codim2PO-sh-mf.jl | 12 ++++++------ examples/codim2PO.jl | 8 ++++---- examples/mittleman.jl | 2 +- examples/pd-1d.jl | 12 ++++++------ src/BifurcationKit.jl | 2 +- src/Continuation.jl | 5 +++-- src/DeflatedContinuation.jl | 4 ++-- src/DeflationOperator.jl | 17 +++++++++-------- src/Newton.jl | 4 ++-- src/NormalForms.jl | 4 ++-- src/bifdiagram/BranchSwitching.jl | 6 +++--- src/codim2/MinAugBT.jl | 2 +- src/codim2/MinAugFold.jl | 2 +- src/codim2/MinAugHopf.jl | 2 +- src/jacobianTypes.jl | 4 ++++ src/periodicorbit/FlowDE.jl | 3 ++- src/periodicorbit/PeriodicOrbitCollocation.jl | 6 ++---- src/periodicorbit/PeriodicOrbitTrapeze.jl | 10 ++++------ src/periodicorbit/PeriodicOrbits.jl | 4 ++-- src/periodicorbit/PoincareRM.jl | 4 ++-- src/periodicorbit/ShootingDE.jl | 4 +++- src/wave/WaveProblem.jl | 2 +- test/simple_continuation.jl | 6 +++--- test/test-cont-non-vector.jl | 12 ++++++------ test/testHopfMA.jl | 6 +++--- test/testJacobianFoldDeflation.jl | 8 ++++---- test/test_newton.jl | 10 +++++----- test/test_wave.jl | 2 +- 37 files changed, 102 insertions(+), 97 deletions(-) diff --git a/examples/SH2d-fronts-cuda.jl b/examples/SH2d-fronts-cuda.jl index cf4a835a..2a2d963a 100644 --- a/examples/SH2d-fronts-cuda.jl +++ b/examples/SH2d-fronts-cuda.jl @@ -118,7 +118,7 @@ prob = BK.BifurcationProblem(F_shfft, AF(sol0), par, (@optic _.l) ; record_from_solution = (x, p; k...) -> norm(x)) opt_new = NewtonPar(verbose = true, tol = 1e-6, linsolver = L, eigsolver = Leig) -sol_hexa = @time newton(prob, opt_new, normN = norminf); +sol_hexa = @time BK.solve(prob, Newton(), opt_new, normN = norminf); println("--> norm(sol) = ", norminf(sol_hexa.u)) plotsol(sol_hexa.u) @@ -126,7 +126,7 @@ plotsol(sol_hexa.u) deflationOp = DeflationOperator(2, 1.0, [sol_hexa.u]) opt_new = @set opt_new.max_iterations = 250 -outdef = @time newton(re_make(prob, u0 = 0.4 .* sol_hexa.u .* AF([exp(-1(x+0lx)^2/25) for x in X, y in Y])), +outdef = @time BK.solve(re_make(prob, u0 = 0.4 .* sol_hexa.u .* AF([exp(-1(x+0lx)^2/25) for x in X, y in Y])), deflationOp, opt_new, normN = x -> maximum(abs, x)) println("--> norm(sol) = ", norm(outdef.u)) plotsol(outdef.u) |> display diff --git a/examples/SH2d-fronts.jl b/examples/SH2d-fronts.jl index 237bf826..32694033 100644 --- a/examples/SH2d-fronts.jl +++ b/examples/SH2d-fronts.jl @@ -59,7 +59,7 @@ optnew = NewtonPar(verbose = true, tol = 1e-8, max_iterations = 20) prob = BifurcationProblem(F_sh, vec(sol0), par, (@optic _.l); J = dF_sh, plot_solution = (x, p; kwargs...) -> (plotsol!((x); label="", kwargs...)),record_from_solution = (x, p; k...) -> (n2 = norm(x), n8 = norm(x, 8)), d2F=d2F_sh, d3F=d3F_sh) # optnew = NewtonPar(verbose = true, tol = 1e-8, max_iterations = 20, eigsolver = EigArpack(0.5, :LM)) -sol_hexa = @time newton(prob, optnew) +sol_hexa = @time BK.solve(prob, Newton(), optnew) println("--> norm(sol) = ", norm(sol_hexa.u, Inf64)) plotsol(sol_hexa.u) @@ -72,7 +72,7 @@ deflationOp = DeflationOperator(2, 1.0, [sol_hexa.u]) optnew = @set optnew.max_iterations = 250 optnewd = @set optnew.max_iterations = 250 -outdef = @time newton( +outdef = @time BK.solve( (@set prob.u0 = 0.4vec(sol_hexa.u) .* vec([exp(-1(x+lx)^2/25) for x in X, y in Y])), deflationOp, # 0.4vec(sol_hexa) .* vec([1 .- exp(-1(x+lx)^2/55) for x in X, y in Y]), optnewd) @@ -129,7 +129,7 @@ end prob2 = @set prob.VF.J = (u, p) -> (du -> dF_sh2(du, u, p)) @reset prob2.u0 = vec(sol0) -sol_hexa = @time newton(prob2, @set optnew.linsolver = ls) +sol_hexa = @time BK.solve(prob2, Newton(), @set optnew.linsolver = ls) println("--> norm(sol) = ", norm(sol_hexa.u, Inf64)) plotsol(sol_hexa.u) ################################################################################################### diff --git a/examples/SH3d.jl b/examples/SH3d.jl index af8801cf..86b15414 100644 --- a/examples/SH3d.jl +++ b/examples/SH3d.jl @@ -123,7 +123,7 @@ prob = BK.BifurcationProblem(F_sh, AF(vec(sol0)), par, (@optic _.l), optnew = NewtonPar(verbose = true, tol = 1e-8, max_iterations = 20, linsolver = @set ls.verbose = 0) @reset optnew.eigsolver = eigSH3d # @reset optnew.linsolver = DefaultLS() -sol_hexa = @time newton(prob, optnew) +sol_hexa = @time BK.solve(prob, Newton(), optnew) println("--> norm(sol) = ", norm(sol_hexa.u, Inf64)) contour3dMakie(sol0) diff --git a/examples/SHpde_snaking.jl b/examples/SHpde_snaking.jl index a91053c6..c207c91c 100644 --- a/examples/SHpde_snaking.jl +++ b/examples/SHpde_snaking.jl @@ -35,7 +35,7 @@ prob = BifurcationProblem(R_SH, sol0, parSH, (@optic _.λ); J = Jac_sp, ) #################################################################################################### optnew = NewtonPar(tol = 1e-12) -sol1 = newton(prob, optnew) +sol1 = BK.solve(prob, Newton(), optnew) Plots.plot(X, sol1.u) opts = ContinuationPar(dsmin = 0.0001, dsmax = 0.01, ds = 0.01, diff --git a/examples/cGL2d-shooting.jl b/examples/cGL2d-shooting.jl index d5babbdc..d953ae80 100644 --- a/examples/cGL2d-shooting.jl +++ b/examples/cGL2d-shooting.jl @@ -118,7 +118,7 @@ prob_sp = SplitODEProblem(f1, f2, sol0_f, (0.0, 120.0), @set par_cgl.r = 1.2; re prob = ODEProblem(Fcgl!, sol0_f, (0.0, 120.0), (@set par_cgl.r = 1.2))#, jac = Jcgl, jac_prototype = Jcgl(sol0_f, par_cgl)) #################################################################################################### # sol = @time solve(prob, Vern9(); abstol=1e-14, reltol=1e-14) -sol = @time solve(prob_sp, ETDRK2(krylov=true); abstol=1e-14, reltol=1e-14, dt = 0.1) #1.78s +sol = @time DifferentialEquations.solve(prob_sp, ETDRK2(krylov=true); abstol=1e-14, reltol=1e-14, dt = 0.1) #1.78s # sol = @time solve(prob, LawsonEuler(krylov=true, m=50); abstol=1e-14, reltol=1e-14, dt = 0.1) # sol = @time solve(prob_sp, CNAB2(linsolve=LinSolveGMRES()); abstol=1e-14, reltol=1e-14, dt = 0.03) diff --git a/examples/cGL2d.jl b/examples/cGL2d.jl index 34e21c55..0780a1b1 100644 --- a/examples/cGL2d.jl +++ b/examples/cGL2d.jl @@ -96,7 +96,7 @@ prob = BK.BifurcationProblem(Fcgl!, vec(sol0), par_cgl, (@optic _.r); J = Jcgl) eigls = EigArpack(1.0, :LM) # eigls = eig_MF_KrylovKit(tol = 1e-8, dim = 60, x₀ = rand(ComplexF64, Nx*Ny), verbose = 1) opt_newton = NewtonPar(tol = 1e-9, verbose = true, eigsolver = eigls, max_iterations = 20) -out = @time newton(prob, opt_newton, normN = norminf) +out = @time solve(prob, Newton(), opt_newton, normN = norminf) #################################################################################################### opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.15, ds = 0.001, p_max = 2.5, detect_bifurcation = 3, nev = 9, plot_every_step = 50, newton_options = (@set opt_newton.verbose = false), max_steps = 1060, n_inversion = 6) br = @time continuation(prob, PALC(), opts_br, verbosity = 0) diff --git a/examples/carrier.jl b/examples/carrier.jl index aa317d36..189c6f6b 100644 --- a/examples/carrier.jl +++ b/examples/carrier.jl @@ -43,7 +43,7 @@ prob = BK.BifurcationProblem(F_carr, zeros(N), par_car, (@optic _.ϵ); record_from_solution) optnew = NewtonPar(tol = 1e-8, verbose = true) -out = @time newton(prob, optnew, normN = norminf) +out = @time BK.solve(prob, Newton(), optnew, normN = norminf) plot(out.u, label = "Solution") optcont = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= -0.01, p_min = 0.05, newton_options = NewtonPar(tol = 1e-8, max_iterations = 20, verbose = false), nev = 40, n_inversion = 6) @@ -67,7 +67,7 @@ function perturbsol(sol, p, id) return sol .+ solp .* sol0 end -outdef1 = @time newton( +outdef1 = @time BK.solve( (@set prob.u0 = perturbsol(-out.u, 0, 0)), deflationOp, # perturbsol(deflationOp[1],0,0), par_def, optdef; diff --git a/examples/chan-af.jl b/examples/chan-af.jl index b1004718..24eca9b6 100644 --- a/examples/chan-af.jl +++ b/examples/chan-af.jl @@ -64,7 +64,7 @@ prob = BifurcationProblem(F_chan, sol0, par_af, (@optic _.α); plot_solution = (x, p; kwargs...) -> plot!(x; label = "l = $(length(x))", kwargs...)) optnew = NewtonPar(tol = 1e-12, verbose = true) -sol = @time BK.newton(prob, optnew, normN = norminf) +sol = @time BK.solve(prob, Newton(), optnew, normN = norminf) plot(sol.u, label="Solution") diff --git a/examples/chan.jl b/examples/chan.jl index b76686f7..11092205 100644 --- a/examples/chan.jl +++ b/examples/chan.jl @@ -27,7 +27,7 @@ prob = BifurcationProblem(F_chan, sol0, par, (@optic _.α); optnewton = NewtonPar(tol = 1e-8, verbose = true) # ca fait dans les 63.59k Allocations -sol = @time newton( prob, optnewton) +sol = @time BK.solve(prob, Newton(), optnewton) optscont = ContinuationPar(dsmin = 0.01, dsmax = 0.5, ds= 0.1, p_max = 4.25, nev = 5, detect_fold = true, plot_every_step = 10, newton_options = NewtonPar(max_iterations = 10, tol = 1e-9, verbose = false), max_steps = 150) @@ -42,8 +42,8 @@ deflationOp = DeflationOperator(2, 1.0, [sol.u]) optdef = NewtonPar(optnewton; tol = 1e-10, max_iterations = 500) -outdef1 = newton(re_make(prob; u0 = sol.u .* (1 .+ 0.01*rand(n))), deflationOp, optdef) -outdef1 = newton(re_make(prob; u0 = sol.u .* (1 .+ 0.01*rand(n))), deflationOp, optdef, Val(:autodiff)) +outdef1 = BK.solve(re_make(prob; u0 = sol.u .* (1 .+ 0.01*rand(n))), deflationOp, optdef) +outdef1 = BK.solve(re_make(prob; u0 = sol.u .* (1 .+ 0.01*rand(n))), deflationOp, optdef, Val(:autodiff)) plot(sol.u, label="newton") plot!(sol0, label="init guess") @@ -51,7 +51,7 @@ plot!(outdef1.u, label="deflation-1") #save newly found point to look for new ones push!(deflationOp, outdef1.u) -outdef2 = @time newton((@set prob.u0 = sol0), deflationOp, optdef; callback = BK.cbMaxNorm(1e5)) +outdef2 = @time BK.solve((@set prob.u0 = sol0), deflationOp, optdef; callback = BK.cbMaxNorm(1e5)) plot!(outdef2.u, label="deflation-2") #################################################################################################### Continuation of the Fold Point using minimally augmented formulation optscont = (@set optscont.newton_options = setproperties(optscont.newton_options; verbose = true, tol = 1e-10)) @@ -110,7 +110,7 @@ P[1,1:2] .= [1, 0.];P[end,end-1:end] .= [0, 1.] lsp = GMRESIterativeSolvers(reltol = 1e-5, N = length(sol.u), restart = 20, maxiter=10, Pl = lu(P)) optnewton_mf = NewtonPar(tol = 1e-9, verbose = false, linsolver = lsp) -out_mf = @time newton(prob2, @set optnewton_mf.verbose = true) +out_mf = @time BK.solve(prob2, Newton(), @set optnewton_mf.verbose = true) plot(brmf, color=:red) @@ -120,7 +120,7 @@ brmf = @time continuation(prob2, PALC(tangent = Bordered(), bls = BorderingBLS(l plot(brmf,color=:blue) alg = brmf.alg -brmf = @time continuation(prob2, MoorePenrose(tangent = alg, directLS = false), (@set opts_cont_mf.newton_options = optnewton_mf)) +brmf = @time continuation(prob2, MoorePenrose(tangent = alg, method = BifurcationKit.iterative), (@set opts_cont_mf.newton_options = optnewton_mf)) plot(brmf,color=:blue) diff --git a/examples/codim2PO-sh-mf.jl b/examples/codim2PO-sh-mf.jl index a03a84c1..f8801982 100644 --- a/examples/codim2PO-sh-mf.jl +++ b/examples/codim2PO-sh-mf.jl @@ -31,9 +31,9 @@ using DifferentialEquations prob_de = ODEProblem(Pop!, z0, (0,600.), par_pop) alg = Rodas5() # alg = Vern9() -sol = solve(prob_de, alg) +sol = DifferentialEquations.solve(prob_de, alg) prob_de = ODEProblem(Pop!, sol.u[end], (0,5.), par_pop, reltol = 1e-10, abstol = 1e-12) -sol = solve(prob_de, Rodas5()) +sol = DifferentialEquations.solve(prob_de, Rodas5()) plot(sol) ################################################################################ @@ -67,7 +67,7 @@ probsh, cish = generate_ci_problem( ShootingProblem(M=3), prob, prob_de, sol, 2. ###### function flow(x0, prob0, tm, p = prob0.p) prob = remake(prob0, u0 = x0, tspan = (0, tm), p = p) - sol = solve(prob, Rodas5()) + sol = DifferentialEquations.solve(prob, Rodas5()) return sol[end] end @@ -97,7 +97,7 @@ AD.pullback_function(AD.FiniteDifferencesBackend(), z -> probsh(z, getparams(pro lspo = GMRESIterativeSolvers(verbose = false, N = length(cish), abstol = 1e-12, reltol = 1e-10) eigpo = EigKrylovKit(x₀ = rand(4)) optnpo = NewtonPar(verbose = true, linsolver = lspo, eigsolver = eigpo) -solpo = newton(probsh, cish, optnpo) +solpo = BK.newton(probsh, cish, optnpo) _sol = BK.get_periodic_orbit(probsh, solpo.u, sol.prob.p) plot(_sol.t, _sol[1:2,:]') @@ -181,8 +181,8 @@ plot!(pd_po_sh, vars = (:ϵ, :b0), branchlabel = "PD") ##### # find the PD NS case par_pop2 = @set par_pop.b0 = 0.45 -sol2 = solve(remake(prob_de, p = par_pop2, u0 = [0.1,0.1,1,0], tspan=(0,1000)), Rodas5()) -sol2 = solve(remake(sol2.prob, tspan = (0,10), u0 = sol2[end]), Rodas5()) +sol2 = DifferentialEquations.solve(remake(prob_de, p = par_pop2, u0 = [0.1,0.1,1,0], tspan=(0,1000)), Rodas5()) +sol2 = DifferentialEquations.solve(remake(sol2.prob, tspan = (0,10), u0 = sol2[end]), Rodas5()) plot(sol2, xlims= (8,10)) probshns, ci = generate_ci_problem( ShootingProblem(M=3), re_make(prob, params = sol2.prob.p), remake(prob_de, p = par_pop2), sol2, 1.; alg = Rodas5(), diff --git a/examples/codim2PO.jl b/examples/codim2PO.jl index 3b47ad35..d3ed704a 100644 --- a/examples/codim2PO.jl +++ b/examples/codim2PO.jl @@ -30,9 +30,9 @@ opts_br = ContinuationPar(p_min = 0., p_max = 20.0, ds = 0.002, dsmax = 0.01, n_ using DifferentialEquations prob_de = ODEProblem(Pop!, z0, (0,600.), par_pop) alg = Rodas5() -sol = solve(prob_de, alg) +sol = DifferentialEquations.solve(prob_de, alg) prob_de = ODEProblem(Pop!, sol.u[end], (0,5.), par_pop, reltol = 1e-8, abstol = 1e-10) -sol = solve(prob_de, Rodas5()) +sol = DifferentialEquations.solve(prob_de, Rodas5()) plot(sol) ################################################################################ @@ -72,7 +72,7 @@ probtrap, ci = generate_ci_problem(PeriodicOrbitTrapProblem(M = 150; jacobian = plot(sol) probtrap(ci, prob.params) |> plot -solpo = newton(probtrap, ci, NewtonPar(verbose = true)) +solpo = BK.newton(probtrap, ci, NewtonPar(verbose = true)) _sol = BK.get_periodic_orbit(probtrap, solpo.u,1) plot(_sol.t, _sol[1:2,:]') @@ -142,7 +142,7 @@ probcoll, ci = generate_ci_problem(PeriodicOrbitOCollProblem(30, 3), prob, sol, plot(sol) probcoll(ci, prob.params) |> plot -solpo = newton(probcoll, ci, NewtonPar(verbose = true)); +solpo = BK.newton(probcoll, ci, NewtonPar(verbose = true)); _sol = BK.get_periodic_orbit(probcoll, solpo.u,1) plot(_sol.t, _sol[1:2,:]') diff --git a/examples/mittleman.jl b/examples/mittleman.jl index ab551140..d2ced859 100644 --- a/examples/mittleman.jl +++ b/examples/mittleman.jl @@ -84,7 +84,7 @@ eigls = EigArpack(20.5, :LM) # eigls = EigKrylovKit(dim = 70) # eigls = EigArpack() opt_newton = NewtonPar(tol = 1e-8, verbose = true, eigsolver = eigls, max_iterations = 20) -sol = newton(prob, opt_newton, normN = norminf) +sol = BK.solve(prob, Newton(), opt_newton, normN = norminf) plotsol(sol.u) #################################################################################################### diff --git a/examples/pd-1d.jl b/examples/pd-1d.jl index 489e584a..b6414df6 100644 --- a/examples/pd-1d.jl +++ b/examples/pd-1d.jl @@ -68,7 +68,7 @@ probBif = BifurcationProblem(Fbr, solc0, par_br, (@optic _.C) ;J = Jbr, eigls = EigArpack(0.5, :LM) optnewton = NewtonPar(eigsolver = eigls, verbose=true, max_iterations = 3200, tol=1e-9) -out = @time newton(probBif, optnewton, normN = norminf) +out = @time BK.solve(probBif, Newton(), optnewton, normN = norminf) plot();plot!(X,out.u[1:N]);plot!(X,solc0[1:N], label = "sol0",line=:dash) optcont = ContinuationPar(dsmax = 0.051, ds = -0.001, p_min = -1.8, detect_bifurcation = 3, nev = 21, plot_every_step = 50, newton_options = optnewton, max_steps = 370, n_inversion = 10, max_bisection_steps = 25) @@ -132,10 +132,10 @@ f1 = DiffEqArrayOperator(par_br.Δ) f2 = NL! prob_sp = SplitODEProblem(f1, f2, solc0, (0.0, 280.0), par_br_hopf) -sol = @time solve(prob_sp, ETDRK2(krylov=true); abstol=1e-14, reltol=1e-14, dt = 0.1, progress = true) +sol = @time DifferentialEquations.solve(prob_sp, ETDRK2(krylov=true); abstol=1e-14, reltol=1e-14, dt = 0.1, progress = true) prob_ode = ODEProblem(Fbr, solc0, (0.0, 280.0), par_br_hopf) -sol = @time solve(prob_ode, Rodas4P(); abstol=1e-14, reltol=1e-7, dt = 0.1, progress = true) +sol = @time DifferentialEquations.solve(prob_ode, Rodas4P(); abstol=1e-14, reltol=1e-7, dt = 0.1, progress = true) orbitsection = Array(sol[:,[end]]) # orbitsection = orbitguess[:, 1] @@ -153,7 +153,7 @@ ls = GMRESIterativeSolvers(reltol = 1e-7, N = length(initpo), maxiter = 50, verb # ls = GMRESKrylovKit(verbose = 0, dim = 200, atol = 1e-9, rtol = 1e-5) optn = NewtonPar(verbose = true, tol = 1e-9, max_iterations = 20, linsolver = ls) # deflationOp = BK.DeflationOperator(2 (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [outpo]) -outposh = @time newton(probSh, initpo, optn; normN = norminf) +outposh = @time BK.newton(probSh, initpo, optn; normN = norminf) BK.converged(outposh) && printstyled(color=:red, "--> T = ", outposh.u[end], ", amplitude = ", BK.getamplitude(probSh, outposh.u, par_br_hopf; ratio = 2),"\n") plot(initpo[1:end-1], label = "Init guess"); plot!(outposh.u[1:end-1], label = "sol") @@ -183,7 +183,7 @@ f2 = NL! prob_sp = SplitODEProblem(f1, f2, solc0, (0.0, 300.0), par_br_pd; abstol=1e-14, reltol=1e-14, dt = 0.01) # solution close to the PD point. -solpd = @time solve(prob_sp, ETDRK2(krylov=true), progress = true) +solpd = @time DifferentialEquations.solve(prob_sp, ETDRK2(krylov=true), progress = true) # heatmap(sol.t, X, sol[1:N,:], color=:viridis, xlim=(20,280.0)) orbitsectionpd = Array(solpd[:,end-100]) @@ -201,7 +201,7 @@ ls = GMRESIterativeSolvers(reltol = 1e-7, N = length(initpo_pd), maxiter = 50, v # ls = GMRESKrylovKit(verbose = 0, dim = 200, atol = 1e-9, rtol = 1e-5) optn = NewtonPar(verbose = true, tol = 1e-9, max_iterations = 12, linsolver = ls) # deflationOp = BK.DeflationOperator(2 (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [outpo]) -outposh_pd = @time newton(BK.set_params_po(probSh,par_br_pd), initpo_pd, optn; +outposh_pd = @time BK.newton(BK.set_params_po(probSh,par_br_pd), initpo_pd, optn; # callback = (state; kwargs...) -> (@show state.x[end];true), normN = norminf) BK.converged(outposh_pd) && printstyled(color=:red, "--> T = ", outposh_pd.u[end], ", amplitude = ", BK.getamplitude(probSh, outposh_pd.u, (@set par_br.C = -0.86); ratio = 2),"\n") diff --git a/src/BifurcationKit.jl b/src/BifurcationKit.jl index 1e5b6327..1b1b8287 100644 --- a/src/BifurcationKit.jl +++ b/src/BifurcationKit.jl @@ -136,7 +136,7 @@ module BifurcationKit export Natural, PALC, Multiple, Secant, Bordered, DefCont, Polynomial, MoorePenrose, MoorePenroseLS, AutoSwitch # newton methods - export NewtonPar, newton, newton_palc, newton_hopf, NonLinearSolution + export NewtonPar, Newton, solve, newton, newton_palc, newton_hopf, NonLinearSolution # continuation methods export ContinuationPar, ContResult, continuation, continuation!, continuation_fold, continuation_hopf, continuation_potrap, eigenvec, eigenvals, get_solx, get_solp, bifurcation_points, SpecialPoint diff --git a/src/Continuation.jl b/src/Continuation.jl index c0eba489..53d62c21 100644 --- a/src/Continuation.jl +++ b/src/Continuation.jl @@ -349,7 +349,7 @@ function Base.iterate(it::ContIterable; _verbosity = it.verbosity) verbose && printstyled("━"^18*" INITIAL GUESS "*"━"^18, bold = true, color = :magenta) # we pass additional kwargs to newton so that it is sent to the newton callback - sol₀ = newton(prob, newton_options; + sol₀ = solve(prob, Newton(), newton_options; normN = it.normC, callback = callback(it), iterationC = 0, @@ -363,7 +363,8 @@ function Base.iterate(it::ContIterable; _verbosity = it.verbosity) verbose && println("──▶ parameter = ", p₀, ", initial step") verbose && printstyled("\n"*"━"^18*" INITIAL TANGENT "*"━"^18, bold = true, color = :magenta) - sol₁ = newton(re_make(prob; params = setparam(it, p₀ + ds / η), u0 = sol₀.u), + sol₁ = solve(re_make(prob; params = setparam(it, p₀ + ds / η), u0 = sol₀.u), + Newton(), newton_options; normN = it.normC, callback = callback(it), diff --git a/src/DeflatedContinuation.jl b/src/DeflatedContinuation.jl index 78b332c3..ffdcca56 100644 --- a/src/DeflatedContinuation.jl +++ b/src/DeflatedContinuation.jl @@ -113,7 +113,7 @@ function updatebranch!(dcIter::DefContIterable, getpredictor!(state, it) pbnew = re_make(it.prob; u0 = _copy(getx(state)), params = setparam(it, current_param)) - sol1 = newton(pbnew, defOp, it.contparams.newton_options, alg.jacobian; + sol1 = solve(pbnew, defOp, it.contparams.newton_options, alg.jacobian; normN = it.normC, callback = it.callback_newton, iterationC = step, @@ -253,7 +253,7 @@ function deflatedContinuation(dcIter::DefContIterable, prob_df = re_make(contIt.prob; u0 = alg.perturb_solution(u0, _p, _idb), params = set(par, lens, _p)) - soln = newton(prob_df, deflationOp, + soln = solve(prob_df, deflationOp, setproperties(optnewton; max_iterations = alg.max_iter_defop); normN = contIt.normC, callback = contIt.callback_newton, diff --git a/src/DeflationOperator.jl b/src/DeflationOperator.jl index bcb94aa8..ba00bec0 100644 --- a/src/DeflationOperator.jl +++ b/src/DeflationOperator.jl @@ -307,7 +307,8 @@ Compared to [`newton`](@ref), the only different arguments are - if passed `Val(:autodiff)`, then `ForwardDiff.jl` is used to compute the jacobian Matrix of the deflated problem - if passed `Val(:fullIterative)`, then a full matrix free method is used for the deflated problem. """ -function newton(prob::AbstractBifurcationProblem, +# GROS PB DE AMBIGUITE LA SEULE SOLUTION EST DE FAIRE solve(prob, Newton(), opt) +function solve(prob::AbstractBifurcationProblem, defOp::DeflationOperator{Tp, Tdot, T, vectype}, options::NewtonPar{T, L, E}, _linsolver::DeflatedProblemCustomLS = DeflatedProblemCustomLS(); @@ -322,25 +323,25 @@ function newton(prob::AbstractBifurcationProblem, # change the linear solver opt_def = @set options.linsolver = linsolver - return newton(deflatedPb, opt_def; kwargs...) + return solve(deflatedPb, Newton(), opt_def; kwargs...) end -function newton(prob::AbstractBifurcationProblem, +function solve(prob::AbstractBifurcationProblem, defOp::DeflationOperator{Tp, Tdot, T, vectype}, options::NewtonPar{T, L, E}, ::Val{:autodiff}; kwargs...) where {Tp, T, Tdot, vectype, L, E} # we create the new functional deflatedPb = DeflatedProblem(prob, defOp, AutoDiff()) - return newton(deflatedPb, options; kwargs...) + return solve(deflatedPb, Newton(), options; kwargs...) end -function newton(prob::AbstractBifurcationProblem, +function solve(prob::AbstractBifurcationProblem, defOp::DeflationOperator{Tp, Tdot, T, vectype}, options::NewtonPar{T, L, E}, ::Val{:fullIterative}; kwargs...) where {Tp, T, Tdot, vectype, L, E} # we create the new functional deflatedPb = DeflatedProblem(prob, defOp, Val(:fullIterative)) - return newton(deflatedPb, options; kwargs...) + return solve(deflatedPb, Newton(), options; kwargs...) end """ @@ -365,11 +366,11 @@ function newton(prob::AbstractBifurcationProblem, linsolver = DeflatedProblemCustomLS(); kwargs...) where {T, vectype, L, E} prob0 = re_make(prob, u0 = x0, params = p0) - sol0 = newton(prob0, options; kwargs...) + sol0 = solve(prob0, Newton(), options; kwargs...) @assert converged(sol0) "Newton did not converge to the trivial solution x0." push!(defOp, sol0.u) prob1 = re_make(prob0, u0 = x1) - sol1 = newton(prob1, defOp, (@set options.max_iterations = 10options.max_iterations), linsolver; kwargs...) + sol1 = solve(prob1, defOp, (@set options.max_iterations = 10options.max_iterations), linsolver; kwargs...) ~converged(sol1) && @error "Deflated Newton did not converge to the non-trivial solution ( i.e. on the bifurcated branch)." @debug "deflated Newton" x0 x1 sol0.u sol1.u # we test if the two solutions are different. We first get the norm diff --git a/src/Newton.jl b/src/Newton.jl index c37393c2..27d3ab35 100644 --- a/src/Newton.jl +++ b/src/Newton.jl @@ -115,7 +115,7 @@ function _newton(prob::AbstractBifurcationProblem, x0, p0, options::NewtonPar; end """ - newton(prob::AbstractBifurcationProblem, options::NewtonPar; normN = norm, callback = (;x, fx, J, residual, step, itlinear, options, x0, residuals; kwargs...) -> true, kwargs...) + solve(prob::AbstractBifurcationProblem, ::Newton, options::NewtonPar; normN = norm, callback = (;x, fx, J, residual, step, itlinear, options, x0, residuals; kwargs...) -> true, kwargs...) This is the Newton-Krylov Solver for `F(x, p0) = 0` with Jacobian w.r.t. `x` written `J(x, p0)` and initial guess `x0`. It is important to set the linear solver `options.linsolver` properly depending on your problem. This linear solver is used to solve ``J(x, p_0)u = -F(x, p_0)`` in the Newton step. You can for example use `linsolver = DefaultLS()` which is the operator backslash: it works well for Sparse / Dense matrices. See [Linear solvers (LS)](@ref) for more informations. @@ -143,7 +143,7 @@ This is the Newton-Krylov Solver for `F(x, p0) = 0` with Jacobian w.r.t. `x` wri !!! warning "Linear solver" Make sure that the linear solver (Matrix-Free...) corresponds to your jacobian (Matrix-Free vs. Matrix based). """ -newton(prob::AbstractBifurcationProblem, options::NewtonPar; kwargs...) = _newton(prob, getu0(prob), getparams(prob), options::NewtonPar; kwargs...) +solve(prob::AbstractBifurcationProblem, ::Newton, options::NewtonPar; kwargs...) = _newton(prob, getu0(prob), getparams(prob), options::NewtonPar; kwargs...) # newton(F, J, x0, p0, options::NewtonPar; kwargs...) = newton(BifurcationProblem(F, x0, p0; J = J), options; kwargs...) # default callback diff --git a/src/NormalForms.jl b/src/NormalForms.jl index 9ab6e36f..4749f4a5 100644 --- a/src/NormalForms.jl +++ b/src/NormalForms.jl @@ -696,7 +696,7 @@ function predictor(bp::NdBranchPoint, δp::T; failures = 0 # we allow for 30 failures of nonlinear deflation while failures < nbfailures - outdef1 = newton(prob, deflationOp, optn, Val(:autodiff); normN = normN) + outdef1 = solve(prob, deflationOp, optn, Val(:autodiff); normN = normN) if converged(outdef1) push!(deflationOp, ampfactor .* outdef1.u) else @@ -757,7 +757,7 @@ function predictor(bp::NdBranchPoint, ::Val{:exhaustive}, δp::T; for ci in igs prob.u0 .= [ci...] * ampfactor # outdef1 = newton(prob, deflationOp, optn, Val(:autodiff); normN = normN) - outdef1 = newton(prob, optn; normN = normN) + outdef1 = solve(prob, Newton(), optn; normN = normN) @debug _ds ci outdef1.converged prob.u0 outdef1.u if converged(outdef1) push!(deflationOp, outdef1.u) diff --git a/src/bifdiagram/BranchSwitching.jl b/src/bifdiagram/BranchSwitching.jl index 78fab432..396c944d 100644 --- a/src/bifdiagram/BranchSwitching.jl +++ b/src/bifdiagram/BranchSwitching.jl @@ -273,9 +273,9 @@ function get_first_points_on_branch(br::AbstractBranchResult, probp = re_make(br.prob; u0 = perturb_guess(bpnf(xsol, ds)), params = setparam(br, bpnf.p + ds)) if usedeflation - solbif = newton(probp, defOpp, optnDf, lsdefop; callback = cbnewton, normN = normn) + solbif = solve(probp, defOpp, optnDf, lsdefop; callback = cbnewton, normN = normn) else - solbif = newton(probp, optnDf; callback = cbnewton, normN = normn) + solbif = solve(probp, Newton(), optnDf; callback = cbnewton, normN = normn) end converged(solbif) && push!(defOpp, solbif.u) end @@ -286,7 +286,7 @@ function get_first_points_on_branch(br::AbstractBranchResult, probm = re_make(br.prob; u0 = perturb_guess(bpnf(xsol, ds)), params = setparam(br, bpnf.p - ds)) if usedeflation - solbif = newton(probm, defOpm, optnDf, lsdefop; callback = cbnewton, normN = normn) + solbif = solve(probm, defOpm, optnDf, lsdefop; callback = cbnewton, normN = normn) else solbif = newton(probm, optnDf; callback = cbnewton, normN = normn) end diff --git a/src/codim2/MinAugBT.jl b/src/codim2/MinAugBT.jl index 642d8848..aceca75d 100644 --- a/src/codim2/MinAugBT.jl +++ b/src/codim2/MinAugBT.jl @@ -332,7 +332,7 @@ function newton_bt(prob::AbstractBifurcationProblem, end # solve the BT equations - sol = newton(prob_bt, optn_bt; normN = normN, kwargs...) + sol = solve(prob_bt, Newton(), optn_bt; normN = normN, kwargs...) # save the solution in BogdanovTakens pbt = get_par_bls(sol.u, 2) diff --git a/src/codim2/MinAugFold.jl b/src/codim2/MinAugFold.jl index ffe31248..588177a0 100644 --- a/src/codim2/MinAugFold.jl +++ b/src/codim2/MinAugFold.jl @@ -229,7 +229,7 @@ function newton_fold(prob::AbstractBifurcationProblem, opt_fold = @set options.linsolver = FoldLinearSolverMinAug() # solve the Fold equations - return newton(prob_f, opt_fold; normN = normN, kwargs...) + return solve(prob_f, Newton(), opt_fold; normN = normN, kwargs...) end function newton_fold(br::AbstractBranchResult, ind_fold::Int; diff --git a/src/codim2/MinAugHopf.jl b/src/codim2/MinAugHopf.jl index f47929c5..61a596ed 100644 --- a/src/codim2/MinAugHopf.jl +++ b/src/codim2/MinAugHopf.jl @@ -237,7 +237,7 @@ function newton_hopf(prob, opt_hopf = @set options.linsolver = HopfLinearSolverMinAug() # solve the hopf equations - return newton(prob_h, opt_hopf, normN = normN, kwargs...) + return solve(prob_h, Newton(), opt_hopf, normN = normN, kwargs...) end function newton_hopf(br::AbstractBranchResult, ind_hopf::Int; diff --git a/src/jacobianTypes.jl b/src/jacobianTypes.jl index f2ec911a..9d715ff8 100644 --- a/src/jacobianTypes.jl +++ b/src/jacobianTypes.jl @@ -1,3 +1,7 @@ +abstract type AbstractNonLinearSolver end + +struct Newton <: AbstractNonLinearSolver end +################################################################################################ abstract type AbstractJacobianType end abstract type AbstractJacobianFree <: AbstractJacobianType end abstract type AbstractJacobianMatrix <: AbstractJacobianType end diff --git a/src/periodicorbit/FlowDE.jl b/src/periodicorbit/FlowDE.jl index 3455cb52..6c09d699 100644 --- a/src/periodicorbit/FlowDE.jl +++ b/src/periodicorbit/FlowDE.jl @@ -1,4 +1,5 @@ -using SciMLBase: remake, solve, ODEProblem, EnsembleProblem, EnsembleThreads, DAEProblem, isinplace as isinplace_sciml +using SciMLBase: remake, ODEProblem, EnsembleProblem, EnsembleThreads, DAEProblem, isinplace as isinplace_sciml +import SciMLBase struct FlowDE{Tprob, Talg, Tjac, TprobMono, TalgMono, Tkwde, Tcb, Tvjp, Tδ} <: AbstractFlow "Store the ODEProblem associated to the flow of the Cauchy problem" diff --git a/src/periodicorbit/PeriodicOrbitCollocation.jl b/src/periodicorbit/PeriodicOrbitCollocation.jl index 8bd6ee46..cfd403e6 100644 --- a/src/periodicorbit/PeriodicOrbitCollocation.jl +++ b/src/periodicorbit/PeriodicOrbitCollocation.jl @@ -933,11 +933,9 @@ function _newton_pocoll(probPO::PeriodicOrbitOCollProblem, prob = WrapPOColl(probPO, jac, orbitguess, getparams(probPO), getlens(probPO), nothing, nothing) if isnothing(defOp) - return newton(prob, options; kwargs...) - # return newton(probPO, jac, orbitguess, par, options; kwargs...) + return solve(prob, Newton(), options; kwargs...) else - # return newton(probPO, jac, orbitguess, par, options, defOp; kwargs...) - return newton(prob, defOp, options; kwargs...) + return solve(prob, defOp, options; kwargs...) end end diff --git a/src/periodicorbit/PeriodicOrbitTrapeze.jl b/src/periodicorbit/PeriodicOrbitTrapeze.jl index b3161943..e664d970 100644 --- a/src/periodicorbit/PeriodicOrbitTrapeze.jl +++ b/src/periodicorbit/PeriodicOrbitTrapeze.jl @@ -861,11 +861,9 @@ function _newton_trap(probPO::PeriodicOrbitTrapProblem, prob = WrapPOTrap(probPO, jac, orbitguess, getparams(probPO.prob_vf), getlens(probPO.prob_vf), nothing, nothing) if isnothing(defOp) - return newton(prob, options; kwargs...) - # return newton(probPO, jac, orbitguess, par, options; kwargs...) + return solve(prob, Newton(), options; kwargs...) else - # return newton(probPO, jac, orbitguess, par, options, defOp; kwargs...) - return newton(prob, defOp, options; kwargs...) + return solve(prob, defOp, options; kwargs...) end else # bordered linear solvers if jacobianPO == :BorderedLU @@ -890,9 +888,9 @@ function _newton_trap(probPO::PeriodicOrbitTrapProblem, prob = WrapPOTrap(probPO, jacPO, orbitguess, getparams(probPO.prob_vf), getlens(probPO.prob_vf), nothing, nothing) if isnothing(defOp) - return newton(prob, (@set options.linsolver = lspo); kwargs...) + return solve(prob, Newton(), (@set options.linsolver = lspo); kwargs...) else - return newton(prob, defOp, (@set options.linsolver = lspo); kwargs...) + return solve(prob, defOp, (@set options.linsolver = lspo); kwargs...) end end end diff --git a/src/periodicorbit/PeriodicOrbits.jl b/src/periodicorbit/PeriodicOrbits.jl index a3c820f9..d5b71583 100644 --- a/src/periodicorbit/PeriodicOrbits.jl +++ b/src/periodicorbit/PeriodicOrbits.jl @@ -212,7 +212,7 @@ function newton(prob::AbstractShootingProblem, kwargs...) jac = _build_jacobian(prob, orbitguess, getparams(prob); δ = δ) probw = WrapPOSh(prob, jac, orbitguess, getparams(prob), lens, nothing, nothing) - return newton(probw, options; kwargs...) + return solve(probw, Newton(), options; kwargs...) end """ @@ -239,7 +239,7 @@ function newton(prob::AbstractShootingProblem, ) where {T, Tp, Tdot, vectype, S, E} jac = _build_jacobian(prob, orbitguess, getparams(prob)) probw = WrapPOSh(prob, jac, orbitguess, getparams(prob), lens, nothing, nothing) - return newton(probw, defOp, options; kwargs...) + return solve(probw, defOp, options; kwargs...) end #################################################################################################### diff --git a/src/periodicorbit/PoincareRM.jl b/src/periodicorbit/PoincareRM.jl index 1145f9e8..185a4363 100644 --- a/src/periodicorbit/PoincareRM.jl +++ b/src/periodicorbit/PoincareRM.jl @@ -119,7 +119,7 @@ function _solve(Π::PoincaréMap{ <: WrapPOSh}, xₛ, par) x₀, par) - solΠ = newton(probΠ, Π.options) + solΠ = solve(probΠ, Newton(), Π.options) ~solΠ.converged && @error "Newton failed!! We did not succeed in computing the Poincaré return map." return solΠ.u end @@ -167,7 +167,7 @@ function _solve(Π::PoincaréMap{ <: WrapPOColl }, xₛ, par) probΠ = BifurcationProblem(mapΠ, x₀, par) - solΠ = newton(probΠ, NewtonPar()) + solΠ = solve(probΠ, Newton(), NewtonPar()) ~solΠ.converged && @error "Newton failed!! We did not succeed in computing the Poincaré return map. Residuals = $(solΠ.residuals)" return solΠ.u end diff --git a/src/periodicorbit/ShootingDE.jl b/src/periodicorbit/ShootingDE.jl index c0db2aff..2d3dba57 100644 --- a/src/periodicorbit/ShootingDE.jl +++ b/src/periodicorbit/ShootingDE.jl @@ -82,7 +82,9 @@ function PoincareShootingProblem(prob::ODEProblem, pSection(out, u, t, integrator) = (hyp(out, u); out .*= integrator.iter > 1) affect!(integrator, idx) = terminate!(integrator) # we put nothing option to have an upcrossing - cb = VectorContinuousCallback(pSection, affect!, hyp.M; interp_points = interp_points, affect_neg! = nothing) + cb = VectorContinuousCallback(pSection, affect!, hyp.M; + interp_points = interp_points, + affect_neg! = nothing) # change ODEProblem -> EnsembleProblem in the parallel case _M = hyp.M parallel = _M == 1 ? false : parallel diff --git a/src/wave/WaveProblem.jl b/src/wave/WaveProblem.jl index 45466e8d..3dec0c67 100644 --- a/src/wave/WaveProblem.jl +++ b/src/wave/WaveProblem.jl @@ -214,7 +214,7 @@ function newton(prob::TWProblem, orbitguess, optn::NewtonPar; kwargs...) jac = (x, p) -> (dx -> prob(x, p, dx)) end probwp = WrapTW(prob, jac, orbitguess, getparams(prob.prob_vf), getlens(prob.prob_vf), record_from_solution(prob.prob_vf), plot_solution(prob.prob_vf)) - return newton(probwp, optn; kwargs...,) + return solve(probwp, Newton(), optn; kwargs...,) end function continuation(prob::TWProblem, diff --git a/test/simple_continuation.jl b/test/simple_continuation.jl index 0146a379..de7cf9b8 100644 --- a/test/simple_continuation.jl +++ b/test/simple_continuation.jl @@ -55,7 +55,7 @@ BK.update(PALC(bls=MatrixBLS(nothing)), ContinuationPar(), nothing).bls.solver = let opts = ContinuationPar(p_min = -3.) prob = BK.BifurcationProblem(F_simple, zeros(10), -1.5, (@optic _); J = Jac_simple) - sol = newton(prob, NewtonPar()) + sol = BK.solve(prob, Newton(), NewtonPar()) iter = ContIterable(prob, PALC(), opts) state = iterate(iter)[1] # jacobian of palc functional @@ -300,8 +300,8 @@ opts = BK.ContinuationPar(dsmax = 0.051, dsmin = 1e-3, ds = 0.001, max_steps = 1 x0 = 0.01 * ones(2) prob = BK.BifurcationProblem(F_simple, x0, -1.5, (@optic _); J = Jac_simple) -x0 = newton(prob, opts.newton_options) -x1 = newton((@set prob.params = -1.45), opts.newton_options) +x0 = BK.solve(prob, Newton(), opts.newton_options) +x1 = BK.solve((@set prob.params = -1.45), Newton(), opts.newton_options) br0 = continuation(prob, PALC(), opts, verbosity=3) BK.get_eigenelements(br0, br0.specialpoint[1]) diff --git a/test/test-cont-non-vector.jl b/test/test-cont-non-vector.jl index fc6b0dfd..3ac64b27 100644 --- a/test/test-cont-non-vector.jl +++ b/test/test-cont-non-vector.jl @@ -26,7 +26,7 @@ prob = BK.BifurcationProblem(F0, [0.8], 1., (@optic _); J = (x, r) -> diagm(0 => 1 .- 3 .* x.^2), Jᵗ = (x, r) -> diagm(0 => 1 .- 3 .* x.^2), d2F = (x, r, v1, v2) -> -6 .* x .* v1 .* v2,) -sol0 = newton(prob, opt_newton0) +sol0 = BK.solve(prob, Newton(), opt_newton0) opts_br0 = ContinuationPar(dsmin = 0.001, dsmax = 0.07, ds= -0.02, p_max = 4.1, p_min = -1., newton_options = setproperties(opt_newton0; max_iterations = 70, tol = 1e-8), detect_bifurcation = 0, max_steps = 150) @@ -73,7 +73,7 @@ sol0 = BorderedArray([0.8], 0.0) opt_newton = NewtonPar(tol = 1e-11, verbose = false, linsolver = linsolveBd()) prob = BK.BifurcationProblem(Fb, sol0, (1., 1.), (@optic _[1]); J = (x, r) -> Jacobian(x, r[1], r[2]), record_from_solution = (x,p;k...) -> x.u[1]) -sol = newton(prob, opt_newton) +sol = BK.solve(prob, Newton(), opt_newton) @test BK.converged(sol) opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= -0.01, p_max = 4.1, p_min = -1., newton_options = setproperties(opt_newton; max_iterations = 70, tol = 1e-8), detect_bifurcation = 0, max_steps = 150, save_sol_every_step = 1) @@ -103,13 +103,13 @@ outfoldco = continuation(br, 1, (@optic _[2]), opts_br; bdlinsolver = BorderingB # test with Newton deflation 1 deflationOp = DeflationOperator(2, 1.0, [zero(sol.u)]) soldef0 = BorderedArray([0.1], 0.0) -soldef1 = newton(BK.re_make(prob, u0 = soldef0), deflationOp, opt_newton) +soldef1 = BK.solve(BK.re_make(prob, u0 = soldef0), deflationOp, opt_newton) push!(deflationOp, soldef1.u) Random.seed!(1231) # test with Newton deflation 2 -soldef2 = newton(BK.re_make(prob, u0 = rmul!(soldef0,rand())), deflationOp, opt_newton) +soldef2 = BK.solve(BK.re_make(prob, u0 = rmul!(soldef0,rand())), deflationOp, opt_newton) #################################################################################################### using KrylovKit @@ -159,7 +159,7 @@ prob = BK.BifurcationProblem(Fr, Jᵗ = (x, p) -> JacobianR(x, p.s), d2F = (x, r, v1, v2) -> RecursiveVec([-6 .* x[ii] .* v1[ii] .* v2[ii] for ii=1:length(x)])) -sol = newton(prob, opt_newton0) +sol = BK.solve(prob, Newton(), opt_newton0) Base.:copyto!(dest::RecursiveVec, in::RecursiveVec) = copy!(dest, in) @@ -194,7 +194,7 @@ outfoldco = continuation(br0sec, 1, (@optic _.s), opts_br0, # try with newtonDeflation # test with Newton deflation 1 deflationOp = DeflationOperator(2, 1.0, [sol.u]) -soldef1 = newton(BK.re_make(prob, u0 = 0.1*(sol.u), params = (r=0., s=1.)), +soldef1 = BK.solve(BK.re_make(prob, u0 = 0.1*(sol.u), params = (r=0., s=1.)), deflationOp, (@set opt_newton0.max_iterations = 20)) push!(deflationOp, soldef1.u) diff --git a/test/testHopfMA.jl b/test/testHopfMA.jl index dfa92f71..bc36c514 100644 --- a/test/testHopfMA.jl +++ b/test/testHopfMA.jl @@ -77,7 +77,7 @@ prob = BifurcationKit.BifurcationProblem(Fbru, sol0, par_bru, (@optic _.l); @test Jbru_sp(sol0, par_bru) - Jbru_ana(sol0, par_bru) |> sparse |> nnz == 0 opt_newton = NewtonPar(tol = 1e-11, verbose = false) -out = newton(BK.re_make(prob, u0 = sol0 .* (1 .+ 0.01rand(2n))), opt_newton) +out = BK.solve(BK.re_make(prob, u0 = sol0 .* (1 .+ 0.01rand(2n))), Newton(), opt_newton) opts_br0 = ContinuationPar(dsmin = 0.001, dsmax = 0.1, ds= 0.01, p_max = 1.8, newton_options = opt_newton, detect_bifurcation = 3, nev = 16, n_inversion = 4) br = continuation(BK.re_make(prob; u0 = out.u, params = (@set par_bru.l = 0.3)), PALC(), opts_br0,) @@ -143,7 +143,7 @@ outhopf = newton(br, 1) pbgopfperso = BK.BifurcationProblem((u, p) -> hopfvariable(u, p), hopfpt, par_bru; J = (x, p) -> Jac_hopf_MA(x, p, hopfvariable),) -outhopf = newton(pbgopfperso, NewtonPar(verbose = false, linsolver = BK.HopfLinearSolverMinAug())) +outhopf = BK.solve(pbgopfperso, Newton(), NewtonPar(verbose = false, linsolver = BK.HopfLinearSolverMinAug())) @test BK.converged(outhopf) # version with analytical Hessian = 2 P(du2) P(du1) QU + 2 PU P(du1) Q(du2) + 2PU P(du2) Q(du1) function d2F(x, p1, du1, du2) @@ -212,7 +212,7 @@ BK.get_time_diff(poTrap, orbitguess_f) # newton to find Periodic orbit _prob = BK.BifurcationProblem((x, p) -> poTrap(x, p), copy(orbitguess_f), (@set par_bru.l = l_hopf + 0.01); J = (x, p) -> poTrap(Val(:JacFullSparse),x,p)) opt_po = NewtonPar(tol = 1e-8, verbose = false, max_iterations = 150) -outpo_f = @time newton(_prob, opt_po) +outpo_f = @time BK.solve(_prob, Newton(), opt_po) @test BK.converged(outpo_f) # println("--> T = ", outpo_f[end]) # flag && printstyled(color=:red, "--> T = ", outpo_f[end], ", amplitude = ", BK.amplitude(outpo_f, n, M; ratio = 2),"\n") diff --git a/test/testJacobianFoldDeflation.jl b/test/testJacobianFoldDeflation.jl index c0bfc9e4..866887bc 100644 --- a/test/testJacobianFoldDeflation.jl +++ b/test/testJacobianFoldDeflation.jl @@ -24,11 +24,11 @@ sol0 = [(i-1)*(n-i)/n^2+0.1 for i=1:n] opt_newton = NewtonPar(tol = 1e-9, verbose = false) prob = BK.BifurcationProblem(F_chan, sol0, (α = 3.3, β = 0.01), (@optic _.α); plot_solution = (x, p; kwargs...) -> (plot!(x;label="", kwargs...))) -out = newton( prob, opt_newton) +out = BK.solve(prob, Newton(), opt_newton) # test with secant continuation opts_br0 = ContinuationPar(dsmin = 0.01, dsmax = 0.15, ds= 0.01, p_max = 4.1, max_steps = 120, newton_options = opt_newton, detect_bifurcation = 3) -br = continuation( prob, PALC(), opts_br0) +br = continuation(prob, PALC(), opts_br0) #################################################################################################### # deflation newton solver, test of jacobian expression deflationOp = DeflationOperator(2, dot, 1.0, [out.u]) @@ -44,9 +44,9 @@ deflationOp = DeflationOperator(2, dot, 1.0, [out.u]) chanDefPb = DeflatedProblem(prob, deflationOp, DefaultLS()) opt_def = NewtonPar(opt_newton; tol = 1e-10, max_iterations = 1000, verbose = false) -outdef1 = newton((@set prob.u0 = out.u .* (1 .+0.01*rand(n))), deflationOp, opt_def) +outdef1 = BK.solve((@set prob.u0 = out.u .* (1 .+0.01*rand(n))), deflationOp, opt_def) # @test BK.converged(outdef1) -outdef1 = newton((@set prob.u0 = out.u .* (1 .+0.01*rand(n))), deflationOp, opt_def, Val(:autodiff)) +outdef1 = BK.solve((@set prob.u0 = out.u .* (1 .+0.01*rand(n))), deflationOp, opt_def, Val(:autodiff)) @test BK.converged(outdef1) #################################################################################################### # Fold continuation, test of Jacobian expression diff --git a/test/test_newton.jl b/test/test_newton.jl index 4c881e8b..6b632f2a 100644 --- a/test/test_newton.jl +++ b/test/test_newton.jl @@ -15,7 +15,7 @@ function test_newton(x0) opts = NewtonPar( tol = Ty(1e-8), verbose = false) prob = BifurcationProblem(F, x0, nothing; J = Jac) - newton(prob, opts; callback = BK.cbMaxNorm(100.0), normN = norminf) + BK.solve(prob, Newton(), opts; callback = BK.cbMaxNorm(100.0), normN = norminf) end #################################################################################################### # we test the regular newton algorithm @@ -140,13 +140,13 @@ end # test of the different newton solvers deflationOp = DeflationOperator(2, dot, _T(1), [[_T(1)]]) -sol = newton(prob, deflationOp, NewtonPar()) +sol = BK.solve(prob, deflationOp, NewtonPar()) @test BK.converged(sol) -sol = newton(prob, deflationOp, NewtonPar(), Val(:autodiff)) +sol = BK.solve(prob, deflationOp, NewtonPar(), Val(:autodiff)) @test BK.converged(sol) -sol = newton(prob, deflationOp, NewtonPar(linsolver = GMRESKrylovKit()),) +sol = BK.solve(prob, deflationOp, NewtonPar(linsolver = GMRESKrylovKit()),) @test BK.converged(sol) -sol = newton(prob, deflationOp, NewtonPar(linsolver = GMRESKrylovKit()), Val(:fullIterative)) +sol = BK.solve(prob, deflationOp, NewtonPar(linsolver = GMRESKrylovKit()), Val(:fullIterative)) @test BK.converged(sol) #################################################################################################### # test newton adapted to branch switching diff --git a/test/test_wave.jl b/test/test_wave.jl index 3c96cba4..87aeae69 100644 --- a/test/test_wave.jl +++ b/test/test_wave.jl @@ -83,7 +83,7 @@ eigls = EigArpack(1.0, :LM) eigls = DefaultEig() # eigls = eig_MF_KrylovKit(tol = 1e-8, dim = 60, x₀ = rand(ComplexF64, Nx*Ny), verbose = 1) opt_newton = NewtonPar(tol = 1e-9, verbose = true, eigsolver = eigls, max_iterations = 20) -out = @time newton(prob, opt_newton, normN = norminf) +out = @time BK.solve(prob, Newton(), opt_newton, normN = norminf) opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.15, ds = 0.001, p_max = 2.5, detect_bifurcation = 3, nev = 9, plot_every_step = 50, newton_options = (@set opt_newton.verbose = false), max_steps = 30, n_inversion = 8, max_bisection_steps=20) br = continuation(prob, PALC(), opts_br, verbosity = 0)