-
-
Notifications
You must be signed in to change notification settings - Fork 33
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
Add BVDAE solvers #217
base: master
Are you sure you want to change the base?
Add BVDAE solvers #217
Conversation
Signed-off-by: Qingyu Qu <2283984853@qq.com>
Signed-off-by: Qingyu Qu <2283984853@qq.com>
Benchmark Results
Benchmark PlotsA plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR. |
Signed-off-by: Qingyu Qu <2283984853@qq.com>
Signed-off-by: Qingyu Qu <2283984853@qq.com>
Signed-off-by: Qingyu Qu <2283984853@qq.com>
Signed-off-by: Qingyu Qu <2283984853@qq.com>
Signed-off-by: Qingyu Qu <2283984853@qq.com>
} | ||
``` | ||
""" | ||
@kwdef struct $(alg){N, J <: BVPJacobianAlgorithm} <: AbstractAscher |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why make each one a different type?
end | ||
|
||
function rkbas!(s, coef, k::Integer, rkb) | ||
# for more guass collocation points per subinterval |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is exactly the radau iia tableau?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The tableaus come from the COLNEW
paper: A New Basis Implementation for a Mixed Order Boundary value ODE solver, which I don't think is the radau iia tableau
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Okay I dug into https://www.researchgate.net/publication/2435083_Collocation_Software_for_Boundary_Value_Differential-Algebraic_Equations a bit and talked with @YingboMa. It's a Gauss-Legendre quadrature https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature. Instad of calling this Ascher we should just call it a FIRK that is GaussLegendre.
Also, I don't really see the reason to be splitting the FIRK methods by different stage numbers, making that a runtime value could mean that we could adapt the order on the fly. But I think we can fix that in the future.
Are there tests with nonlinearleastsquaresproblem? I didn't see any. |
Haven't found an appropriate test example for its case, previous work on this PR only ensures the "standard" BVDAE part is working fine, but I suspect COLDAE support over-constraints or under-constraints BVDAE |
Trying to understand the general system here. I think it's fine for the Ascher (Gauss-Legendre) to be kept separate from the other FIRK because they have a lot going on with their adaptivity algorithm and such, that's fine. Maybe we keep the name Ascher because of that but in the docstring mention it's Gauss-Legendre but with Ascher's adaptivity algorithm. But the bigger point, for BVDAEs, is there anything that limits them to the GL methods? What I ask is, any FIRK method should work for the DAE in mass matrix form right? It would be good to test that this is method independent. I'm trying to see where this mass matrix is being used and I only see https://github.com/SciML/BoundaryValueDiffEq.jl/pull/217/files#diff-f206170122c6dd64b293fa8ee49c4b402d313b6ad26dfccc7bb1ad11624f91f9R65 so I'm a bit confused. |
Good points, will add them then.
Yeah, I can test if the FIRK tableaus can plug into this implementation but I suspected there would be extra work comparing that we can just replace FIRK tableaus with MIRK to get FIRK collocation methods for common BVPs. The current implementation follows the pattern of COLDAE, in which IIUC only needs to specify the differential vars, so the mass matrix here is only used for stating which equations should be the differential or algebraic equations. |
Signed-off-by: Qingyu Qu <2283984853@qq.com>
The BVDAE solvers are almost done now, but need further optimization in the style and improving its performance. The BVDAE solvers are similar to MIRK and FIRK methods, and their numerical workflow can be roughly divided into several steps:
(1) Build collocation equations with RK tableau. With Runge-Kutta coefficients at collocation points, build the collocation equations at each subinterval. After the collocation equations are built, we need to use AlmostBlockDiagonals.jl) to solve the built linear systems, which has an almost block diagonal form, we can just factorize(pLU) and substitute to find the solutions.
(2) Construct a nonlinear system. Construct a nonlinear system using collocation equations. In the nonlinear solving process, residual control is deployed to obtain the numerical solution. The sparsity of the nonlinear system is also configured according to what we did in MIRK and FIRK.
(3) Error estimation. While we use defect control for error control in MIRK and FIRK, BVDAE solvers utilize a common error control mechanism, which compares the two solutions obtained from different mesh(solutions from original mesh and halved mesh) and estimates the errors.
(4) Refine mesh according to the error norm and degree of equidistribution. With the error estimation we get, we compute the degree of equidistribution of the approximate integral of function and accordingly decide to refine the mesh or just halve the mesh.
(5) Proceed for next iteration or successful return.
There are still several issues that need to be fixed:
BVDAEProblem
like what we did forDAEProblem
Also implemented an interface for COLDAE in ODEInterface.jl, since the ODEInterface.jl is not actively maintained, just put the implementation in fork: https://github.com/ErikQQY/ODEInterface.jl
Fix: #76