diff --git a/src/PALC/continuation.jl b/src/PALC/continuation.jl index 92819ce..ac0de9f 100644 --- a/src/PALC/continuation.jl +++ b/src/PALC/continuation.jl @@ -25,7 +25,7 @@ function continuation( solvers = PALCSolverCache(p, alg, cache, newton_iter, newton_tol, newton_max_resid) # Initialize continuation - initialize_palc!(initial_tangent, cache, alg, p, solvers, trace) + initialize_palc!(initial_tangent, cache, alg, p, solvers, trace, dsmin, dsmax) # Continuation loop continuation!( diff --git a/src/PALC/initialization.jl b/src/PALC/initialization.jl index a1b3f2c..e2a4876 100644 --- a/src/PALC/initialization.jl +++ b/src/PALC/initialization.jl @@ -46,38 +46,46 @@ function compute_initial_tangent( prob::ContinuationProblem, solvers, trace::AbstractTraceLevel, -) - # Perturb the natural continuation parameter to construct initial tangent - λ_pert = method.step_scale_factor * cache.ds - perturb_natural_continuation_parameter!(cache, λ_pert) + dsmin, + dsmax +) + done = false + while !done + # Perturb the natural continuation parameter to construct initial tangent + λ_pert = method.step_scale_factor * cache.ds + perturb_natural_continuation_parameter!(cache, λ_pert) - # Print trace if desired - print_initialization_trace(cache, trace, 2) + # Print trace if desired + print_initialization_trace(cache, trace, 2) - # Resolve the problem with perturbed parameter - u_sol, retcode = solve_natural_nlp!(solvers, cache.u0, trace) + # Resolve the problem with perturbed parameter + u_sol, retcode = solve_natural_nlp!(solvers, cache.u0, trace) - if SciMLBase.successful_retcode(retcode) - # Compute and set the initial tangent - u0 = cache.u0 # The current solution (from first solve) - n = length(u0) - δuλ0 = cache.δuλ0 # The predicted tangent direction - - # Set secant direction - δuλ0[1:n] .= u_sol .- cache.u0 - δuλ0[end] = λ_pert - - # Normalize the tangent direction and ensure it is in the positive λ direction - sf = λ_pert < 0 ? -1.0 : 1.0 - n_inv = sf / norm(δuλ0) - δuλ0 .*= n_inv - - update_tangent!(cache, δuλ0) # Update the tangent direction in cache - cache.δuλ0_initial .= δuλ0 # Save initial tangent - else - error( - "Solve to compute initial tangent failed! Consider reducing perturbation size." - ) + if SciMLBase.successful_retcode(retcode) + # Compute and set the initial tangent + u0 = cache.u0 # The current solution (from first solve) + n = length(u0) + δuλ0 = cache.δuλ0 # The predicted tangent direction + + # Set secant direction + δuλ0[1:n] .= u_sol .- cache.u0 + δuλ0[end] = λ_pert + + # Normalize the tangent direction and ensure it is in the positive λ direction + sf = λ_pert < 0 ? -1.0 : 1.0 + n_inv = sf / norm(δuλ0) + δuλ0 .*= n_inv + + update_tangent!(cache, δuλ0) # Update the tangent direction in cache + cache.δuλ0_initial .= δuλ0 # Save initial tangent + done = true + elseif cache.ds == dsmin # we have tried with the smallest allowable ds and still failed, so error. + error( + "Solve to compute initial tangent failed! Consider reducing dsmin." + ) + else + scale_and_clamp_ds!(cache, 0.5, dsmin, dsmax) # we are not yet at ds=dsmin, so scale and re-try + end end return nothing end @@ -89,6 +97,8 @@ function compute_initial_tangent( prob::ContinuationProblem, solvers, trace::AbstractTraceLevel, + dsmin, + dsmax ) # Get user tangent n = length(cache.u0) @@ -114,6 +124,8 @@ function compute_initial_tangent( prob::ContinuationProblem, solvers, trace::AbstractTraceLevel, + dsmin, + dsmax ) n = length(cache.u0) uλpred = cache.uλpred @@ -179,6 +191,8 @@ function initialize_palc!( p::ContinuationProblem, solvers, trace::AbstractTraceLevel, + dsmin, + dsmax ) # Get information from the problem u0 = p.u0 @@ -206,7 +220,7 @@ function initialize_palc!( end # Compute the initial tangent - compute_initial_tangent(initial_tangent, cache, alg, p, solvers, trace) + compute_initial_tangent(initial_tangent, cache, alg, p, solvers, trace, dsmin, dsmax) return nothing end