Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Pseudo-Transient Method #228

Closed
wants to merge 8 commits into from
Closed

[WIP] Pseudo-Transient Method #228

wants to merge 8 commits into from

Conversation

yonatanwesen
Copy link
Contributor

@yonatanwesen yonatanwesen commented Oct 4, 2023

This PR adds pseudo-transient methods using adaptive step size described here created this new PR since #215 had too many merge conflicts

Comment on lines 409 to 460
@testset "PseudoTransient" begin

#iip
function test_f!(du, u, p)
du[1] = 2 - 2u[1]
du[2] = u[1] - 4u[2]
return du
end

#oop
simple_test(u, p) = 0.9f0 .* u .- u

#test jacobian free PseudoTransient method
function f!(res, u, p)
@. res = u * u - p
end

@testset "[IIP] u0: $(typeof(u0))" for u0 in (zeros(2),)
probN = NonlinearProblem{true}(test_f!, u0)
sol = solve(probN, PseudoTransient())
@test sol.retcode == ReturnCode.Success
du = zeros(2)
test_f!(du, sol.u, 4.0)
@test du≈[0, 0] atol=1e-7
end

@testset "[OOP] u0: $(typeof(u0))" for u0 in ([1.0f0],)
probN = NonlinearProblem{false}(simple_test, u0)
sol = solve(probN, PseudoTransient(alpha_initial = 1.0))
@test sol.retcode == ReturnCode.Success
@test abs(sol.u[1]) <= 1.0f-4
end

@testset "Jacobian Free PT" begin
u0 = [1.0, 1.0]
p = 2.0
probN = NonlinearProblem(f!, u0, p)
linsolve = LinearSolve.KrylovJL_GMRES()
sol = solve(probN,
PseudoTransient(alpha_initial = 10.0, linsolve = linsolve),
reltol = 1e-9)
@test sol.retcode == ReturnCode.Success
end

@testset "NewtonRaphson Fails" begin
p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0]
probN = NonlinearProblem{false}(newton_fails, u0, p)
sol = solve(probN, PseudoTransient(alpha_initial = 1.0), abstol = 1e-10)
@test all(abs.(newton_fails(sol.u, p)) .< 1e-10)
end
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use similar setup as the NewtonRaphson tests.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants