WIP: Implement the ROCK2 method #374
Conversation
2 similar comments
src/tableaus/rkc_tableaus.jl
Outdated
| ms::SVector{46, Int} | ||
| fp1::SVector{46, T} | ||
| fp2::SVector{46, T} | ||
| recfi::SVector{4476, T2} |
There was a problem hiding this comment.
This should be a Vector. You typically don't want to use SVector for large arrays. All of these should probably be arrays, but some testing would need to be done on the 46.
src/tableaus/rkc_tableaus.jl
Outdated
| recfi::SVector{4476, T2} | ||
| end | ||
|
|
||
| Base.@pure function ROCK2ConstantCache{T,T2}(::Type{T},::Type{T2}) |
There was a problem hiding this comment.
What code generated these coefficients?
|
Looks like a great start! Could you link what you're getting this from? It looks like there's a lot of hardcoded constants, so if this comes from another code we need to be careful about the license. It should be an MIT-compatible license. If it's generated from some formula in the paper, the generated code should be included (commented) and cited. |
|
The hardcoded constants are from http://www.unige.ch/~hairer/software.html |
|
Cool, that page is all MIT licensed: http://www.unige.ch/~hairer/prog/licence.txt so we're good. |
Codecov Report
@@ Coverage Diff @@
## master #374 +/- ##
==========================================
- Coverage 95.93% 94.98% -0.95%
==========================================
Files 74 78 +4
Lines 9073 8933 -140
==========================================
- Hits 8704 8485 -219
- Misses 369 448 +79
Continue to review full report at Codecov.
|
src/perform_step/rkc_perform_step.jl
Outdated
| temp2 = 1 - temp3 | ||
| ci1 = temp1 + temp2 * ci2 + temp3 * ci3 | ||
| u = temp1 * u + temp2 * gprev + temp3 * gprev2 | ||
| i < ms[cache.mdeg] && (gprev2 = gprev, gprev = u) |
There was a problem hiding this comment.
I think (gprev2 = gprev, gprev = u) should be (gprev2 = gprev; gprev = u).
There was a problem hiding this comment.
Also, the conditioning in here seems superfluous. Nevermind.
src/perform_step/rkc_perform_step.jl
Outdated
| temp3 = temp2 * ( u - gprev2) | ||
| u = gprev + temp1 * u + temp3 | ||
| ci1 = max(abs(u), abs(uprev)) * rto | ||
| calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm) |
There was a problem hiding this comment.
Use calculate_residuals for the mutable version.
src/perform_step/rkc_perform_step.jl
Outdated
| if ntol = 0 # what is ntol? | ||
| temp3 = temp2 * ( u - gprev2) | ||
| u = gprev + temp1 * u + temp3 | ||
| ci1 = max(abs(u), abs(uprev)) * rto |
There was a problem hiding this comment.
This doesn't look right. u and uprev are vectors in here.
src/perform_step/rkc_perform_step.jl
Outdated
| gprev = u + temp1 * gprev2 | ||
| ci1 += temp1 | ||
| # error estimate | ||
| if ntol = 0 # what is ntol? |
There was a problem hiding this comment.
src/perform_step/rkc_perform_step.jl
Outdated
| ms[cache.mdeg] < 2 && u = gprev | ||
| # for the second to the ms[cache.mdeg] th stages | ||
| for i in 2:ms[cache.mdeg] | ||
| temp1 = h * recf[cache.recind + 2 * (i - 2) + 1] |
src/perform_step/rkc_perform_step.jl
Outdated
| temp2 = 1 - temp3 | ||
| ci1 = temp1 + temp2 * ci2 + temp3 * ci3 | ||
| u = temp1 * u + temp2 * gprev + temp3 * gprev2 | ||
| i < ms[cache.mdeg] && (gprev2 = gprev, gprev = u) |
There was a problem hiding this comment.
Also, the conditioning in here seems superfluous. Nevermind.
src/perform_step/rkc_perform_step.jl
Outdated
| ci3 = t | ||
| gprev2 = uprev | ||
| gprev = uprev + temp1 * fsalfirst | ||
| ms[cache.mdeg] < 2 && u = gprev |
src/perform_step/rkc_perform_step.jl
Outdated
| mdeg = Int(floor(sqrt((1.5 + dt * integrator.eigen_est)/0.811) + 1)) | ||
| if mdeg >= 200 | ||
| h = 0.8 * (200 ^ 2 * 0.811 - 1.5)/integrator.eigen_est | ||
| dt = 0.8 * (200 ^ 2 * 0.811 - 1.5)/integrator.eigen_est |
There was a problem hiding this comment.
this won't actually update the dt. A more sophisticated approach will be required here. Is the dt only ever decreased by this?
|
How's this doing? |
|
Looks like there's been a lot of progress and a lot of it looks great. However, I don't understand the implementation of your core inner loop there. What paper did you work from? Looking at http://www.dumkaland.org/publications/content.pdf , I don't understand what the mu, nu, and kappa terms are, or how your loop corresponds to (26) and (27). |
YingboMa
left a comment
YingboMa
left a comment
There was a problem hiding this comment.
We need more tests to ensure the correctness of this algorithm, as there are many recurrence parameters and a simple linear convergence test only covers a few of them.
|
So the Then is I think I am getting it but this can be cleaned up a bit so that way it's clearer at first glance what's going on. First off, why not split this up into Other than that, I'm not entirely sure how you're choosing the maximum stages. You have the eigenvalue computation here https://github.com/tyfff/OrdinaryDiffEq.jl/blob/ebff1738d6dfe08d7c6277d80abda1cde842ce3e/src/rkc_utils.jl#L4 but I don't see it actually used? Could you describe the stage number choice? |
|
That test failure is fixed on master |
|
This is a pretty straight port of the Fortran code and it passes tests so I decided to merge. I still want to do a reorganization so that way it's much easier to follow, but this will be a separate issue and PR. Thanks! This is a very strong first contribution! |
4 participants
I am new to Julia, programming, and open-source. If there is any mistake, please tell me. @YingboMa is currently mentoring me.