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

GeneralDomain documentation is unclear #133

Open
DanielVandH opened this issue Jan 14, 2023 · 2 comments
Open

GeneralDomain documentation is unclear #133

DanielVandH opened this issue Jan 14, 2023 · 2 comments

Comments

@DanielVandH
Copy link
Member

DanielVandH commented Jan 14, 2023

The documentation for GeneralDomain shows the docstring

GeneralDomain(g, u=nothing; nlsolve=NLSOLVEJL_SETUP(), save=true,
                       abstol=nothing, scalefactor=nothing, autonomous=numargs(g)==2,
                       nlopts=Dict(:ftol => 10*eps()))

which should now be

GeneralDomain(g, u = nothing; nlsolve = NLSOLVEJL_SETUP(), save = true,
                       abstol = nothing, scalefactor = nothing,
                       autonomous = maximum(SciMLBase.numargs(g)) == 3,
                       nlopts = Dict(:ftol => 10 * eps()))

Moreover, it initially says that g should be of the form g(u, resid) or g(t, u, resid), and then later says it should be g(resid, u, p) or g(resid, u, p, t). What should it be? It's not exactly clear to me what the "residuals" are either: a vector that tests each component of the solution? A single value testing all components? I thought it would be the former, but the documentation says

- `g`: ... which is zero when the value is in the domain.

which suggests the latter. I guess judging from

function isaccepted(u, p, t, abstol, f::GeneralDomainAffect{autonomous, F, T, S, uType},
resid) where {autonomous, F, T, S, uType}
# calculate residuals
if autonomous
f.g(resid, u, p)
else
f.g(resid, u, p, t)
end
# accept time step if residuals are smaller than the tolerance
if typeof(abstol) <: Number
all(x -> x < abstol, resid)
else
# element-wise comparison
length(resid) == length(abstol) ||
throw(DimensionMismatch("numbers of residuals and tolerances do not match"))
all(x < y for (x, y) in zip(resid, abstol))
end
end

it should be that g is g(resid, u, p) or g(resid, u, p, t), and resid could actually be any mutable value, not necessarily with length(resid) == length(u) (doesn't seem like resid can be defined manually directly actually, and typeof(resid) must == typeof(u)).

I'm also guessing that p is just integrator.p, or is there a way to add on extra parameters here without a closure?

@ChrisRackauckas
Copy link
Member

Moreover, it initially says that g should be of the form g(u, resid) or g(t, u, resid), and then later says it should be g(resid, u, p) or g(resid, u, p, t). What should it be?

The latter.

It's not exactly clear to me what the "residuals" are either: a vector that tests each component of the solution?

Yes, the residuals of the nonlinear system.

I'm also guessing that p is just integrator.p, or is there a way to add on extra parameters here without a closure?

That's all it is.

Though, we should probably extend this a little bit and change the nonlinear solving to use NonlinearSolve.jl. That would be a breaking change, but quite worth it because this code is really old. Make it so that way f always has to have t to simplify the code too.

ChrisRackauckas added a commit that referenced this issue Jan 22, 2023
Clarify GeneralDomain docstring (#133)
@ChrisRackauckas
Copy link
Member

@ArnoStrouwen can you tackle this one?

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

No branches or pull requests

2 participants