Add defect control adaptivity #89

ErikQQY · 2023-05-29T13:04:04Z

Fix: #12
Fix: #54

Relating reference: https://netlib.org/ode/mirkdc.f

I just built the basic structure for the defect control adaptivity, it still needs a lot of work to get this program working.

The detailed explanation can be viewed on my blog post

ping @ChrisRackauckas @YingboMa @avik-pal

Signed-off-by: ErikQQY <2283984853@qq.com>

src/solve.jl

Signed-off-by: ErikQQY <2283984853@qq.com>

…Eq.jl

Signed-off-by: ErikQQY <2283984853@qq.com>

ChrisRackauckas · 2023-06-17T16:05:06Z

src/solve.jl

+function redistribute(mesh_current::Vector,
+                      n::Int64,
+                      Nsub_star::Int64,
+                      s_hat::Vector{Float64})
+    mesh_new = zeros(Float64, Nsub_star + 1)
+    sum = 0.0
+    for k in 1:n
+        sum += s_hat[k] * (mesh_current[k + 1] - mesh_current[k])
+    end
+    zeta = sum / Nsub_star
+    k::Int64 = 1
+    i::Int64 = 0
+    mesh_new[1] = mesh_current[1]
+    t = mesh_current[1]
+    integral::Float64 = 0.0
+    while k <= n
+        next_piece = s_hat[k] * (mesh_current[k + 1] - t)
+        if (integral + next_piece) > zeta
+            mesh_new[i + 2] = (zeta - integral) / s_hat[k] + t
+            t = mesh_new[i + 2]
+            i += 1
+            integral = 0
+        else
+            integral += next_piece
+            t = mesh_current[k + 1]
+            k += 1
+        end
+    end
+    mesh_new[Nsub_star + 1] = mesh_current[n + 1]
+    return mesh_new
+end


Can you describe what this one is doing?

The redistribute part generates a new mesh to make sure that on the next iteration, the new solution on the new mesh would have a smaller defect estimation.

Based on the judging quantities we computed before using defect and mesh function values, we compare the ratio $\frac{r1}{r3}$ with $\rho$ to decide whether we halve the mesh or make a redistribute.

s_hat = zeros(Float64, n) for i in 1:n h = mesh_current[i + 1] - mesh_current[i] norm = abs(defect[i, idamax(defect[i, :])]) s_hat[i] = (norm / tol)^(1.0 / (p + 1)) / h if s_hat[i] * h > r1 r1 = s_hat[i] * h end r2 = r2 + s_hat[i] * h end r3 = r2 / n n_predict::Int64 = round(Int, (safety_factor * r2) + 1) if abs((n_predict - n) / n) < 0.1 n_predict = round(Int, 1.1 * n) end

If we need to redistribute(r1/r3>rho), we redistribute a new mesh, I keep some detailed notes on the redistribute process here: https://erikqqy.github.io/2023/06/04/gsoc-1st-week/#Redistribute-mesh

ChrisRackauckas · 2023-06-17T16:05:53Z

src/solve.jl

+    safety_factor = 1.3
+    rho = 1.0 # Set rho=1 means mesh distribution will take place everytime.
+    upper_new_mesh = 4.0
+    lower_new_mesh = 0.5
+    r1 = 0.0
+    r2 = 0.0
+    Nsub_star = 0
+    info = 0


are these pulled from the paper?

The BVP book does have the specification of r1, r2 and r3

But both the book and the paper don't elaborate on the setting of the other parameters, so I set these parameters according to the FORTRAN routine.

ChrisRackauckas · 2023-06-17T16:06:42Z

src/solve.jl

+            info = -1
+        else
+            println("Mesh redistributing")
+            mesh_new = redistribute(mesh_current, n, Nsub_star, s_hat)


I don't undersatnd why one would redistribute. Wouldn't this make it so you aren't clustering the points where the error is the most anymore?

From my understanding, after we redistribute the mesh and by using the new mesh to generate a new discrete solution, the new solution would have a smaller error.

ChrisRackauckas · 2023-06-17T16:07:03Z

src/solve.jl

+
+        # Change the original solution Vector{Vector{Float64}} to a matrix
+        # has size of (number of subintervals, length of equation)
+        global Y = transpose(reduce(hcat, S.y))


why is this global?

The discrete solution Y is a local variable in the loop of iteration, so I set Y as global so we can finally build a solution based on Y.

Signed-off-by: ErikQQY <2283984853@qq.com>

ChrisRackauckas · 2023-06-24T16:15:32Z

What's left to complete this? Looks like there's some merge issues that are just from the formatter.

ErikQQY · 2023-06-24T16:54:39Z

The main issue here is that the defect norm from the defect estimation part is always larger than the threshold, making the redistribution part hit an endless loop.

I think that the main problem can be located in the non-linear solving part of BundaryValueDiffEq.jl and MIRKDC, they are using different Newton iteration techniques, which leads to the difference in discrete stages k_discrete, thus the defect we finally get maybe result in a different outcome.

Signed-off-by: ErikQQY <2283984853@qq.com>

avik-pal · 2023-07-10T15:15:27Z

src/cache.jl

@@ -17,22 +17,26 @@ function alg_cache{T,U}(alg::MIRK4, S::BVPSystem{T,U})
 end
 =#

-struct MIRK4GeneralCache{kType} <: GeneralMIRKCache


Do we need this to be mutable? It seems that the array elements are being updated later.

for i in 1:(N - 1) h = x[i + 1] - x[i] # Update K update_K!(S, cache, TU, i, h) for r in 1:s x_new = x[i] + c[r] * h y_new = (1 - v[r]) * y[i] + v[r] * y[i + 1] if r > 1 y_new += h * sum(j -> X[r, j] * K[j], 1:(r - 1)) end temp = zeros(Float64, M) fun!(temp, y_new, S.p, x_new) K[r] = temp[:] end # Update residual residual[i] = y[i + 1] - y[i] - h * sum(j -> b[j] * K[j], 1:s) cache.k_discrete[i, :] = K[:] end

The MIRK cache stores the discrete stages from evaluating PHI in the nonlinear solving process. Since we only need discrete stages from the final nonlinear solving step, I think the previous steps can be omitted, but I didn't figure out an executable way, e.g. judging if this step is the last time we evaluate PHI🤔.

It's mutating the values but not the cache itself

Signed-off-by: ErikQQY <2283984853@qq.com>

ChrisRackauckas · 2023-07-18T14:59:06Z

src/adaptivity.jl

+        return length(mesh) - 1
+    else
+        ind = findfirst(x -> x >= t, mesh)
+        i::Int64 = copy(ind)


why does this need a declaration?

oh the nothing case. Ehh that's kind of not nice but it works

ChrisRackauckas · 2023-07-18T15:00:16Z

src/adaptivity.jl

+    MxNsub = 3000
+
+    safety_factor = 1.3
+    rho = 1.0 # Set rho=1 means mesh distribution will take place everytime.
+    upper_new_mesh = 4.0
+    lower_new_mesh = 0.5
+    r1 = 0.0
+    r2 = 0.0
+    Nsub_star = 0


Surface things as kwargs in a follow up?

src/adaptivity.jl

ChrisRackauckas · 2023-07-18T15:01:17Z

src/adaptivity.jl

+            # Mesh redistribution fails
+            #println("New mesh would be too large")
+            info = ReturnCode.Failure
+            mesh_new = Nothing


Don't use the datatype, that won't specialize

ChrisRackauckas · 2023-07-18T15:01:45Z

src/adaptivity.jl

+Generate a new mesh.
+"""
+function redistribute(mesh_current::Vector,
+    n::Int64,


Why require Int64?

ChrisRackauckas · 2023-07-18T15:02:23Z

src/adaptivity.jl

+            k += 1
+        end
+    end
+    mesh_new[Nsub_star + 1] = mesh_current[n + 1]


This will fail on 32-bit OS with the coercion to Int64.

src/adaptivity.jl

ChrisRackauckas · 2023-07-18T15:03:04Z

src/adaptivity.jl

+    mesh_new[1] = mesh_current[1]
+    for i in 1:n
+        mesh_new[2 * i + 1] = mesh_current[i + 1]
+        mesh_new[2 * i] = (mesh_current[i + 1] + mesh_current[i]) / 2.0


Suggested change

mesh_new[2 * i] = (mesh_current[i + 1] + mesh_current[i]) / 2.0

mesh_new[2 * i] = (mesh_current[i + 1] + mesh_current[i]) / 2

src/adaptivity.jl

ChrisRackauckas · 2023-07-18T15:04:12Z

src/adaptivity.jl

+    f_sample_1, f_sample_2 = zeros(Float64, len), zeros(Float64, len)
+    def_1, def_2 = zeros(Float64, len), zeros(Float64, len)
+    temp_1, temp_2 = zeros(Float64, len), zeros(Float64, len)
+    estimate_1, estimate_2 = zeros(Float64), zeros(Float64)


Assumes Float64

ChrisRackauckas · 2023-07-18T15:04:28Z

src/adaptivity.jl

+    # Evaluate at the first sample point
+    weights_1, weights_1_prime = interp_weights(tau_star, alg)
+    # Evaluate at the second sample point
+    weights_2, weights_2_prime = interp_weights(1.0 - tau_star, alg)


Suggested change

weights_2, weights_2_prime = interp_weights(1.0 - tau_star, alg)

weights_2, weights_2_prime = interp_weights(1 - tau_star, alg)

etc. shouldn't assume float64

I think while we can know the type of tau_star, why don't we directly use Float to avoid redundant type conversion?

ChrisRackauckas · 2023-07-18T15:04:49Z

src/adaptivity.jl

+function interp_weights(tau, alg::Union{GeneralMIRK, MIRK})
+    if alg_order(alg) == 4
+        t2 = tau * tau
+        tm1 = tau - 1.0


Suggested change

tm1 = tau - 1.0

tm1 = tau - 1

Etc., shouldn't assume Float64

src/adaptivity.jl

Signed-off-by: ErikQQY <2283984853@qq.com>

ChrisRackauckas · 2023-07-22T21:18:53Z

Needs to rebase for the MIRK5, then I think we can merge and continue some of this in a follow up.

Signed-off-by: ErikQQY <2283984853@qq.com>

ErikQQY · 2023-07-24T15:33:17Z

@ChrisRackauckas bump

ChrisRackauckas · 2023-07-24T15:36:20Z

Open issues for the points that weren't addressed, but let's go for it. 🎉 🎉 🎉 🎉 🎉 🎉 🎉 🎉 🎉

Add defect control adaptivity

e2fc944

Signed-off-by: ErikQQY <2283984853@qq.com>

ErikQQY marked this pull request as draft May 29, 2023 13:04

avik-pal requested changes May 30, 2023

View reviewed changes

src/solve.jl Outdated Show resolved Hide resolved

ErikQQY added 8 commits June 1, 2023 00:17

Fix comments

33b221e

Signed-off-by: ErikQQY <2283984853@qq.com>

format

f2dd8f6

Signed-off-by: ErikQQY <2283984853@qq.com>

Fix errors

7d028ba

Signed-off-by: ErikQQY <2283984853@qq.com>

Fix more errors

adcc35a

Signed-off-by: ErikQQY <2283984853@qq.com>

Make mesh refinement work

8863e45

Signed-off-by: ErikQQY <2283984853@qq.com>

format

9f9e437

Signed-off-by: ErikQQY <2283984853@qq.com>

Ready for review

1bfc50e

Signed-off-by: ErikQQY <2283984853@qq.com>

Merge branch 'master' of https://github.com/ErikQQY/BoundaryValueDiff…

642e5c8

…Eq.jl

ErikQQY marked this pull request as ready for review June 11, 2023 08:09

ErikQQY and others added 3 commits June 11, 2023 16:21

format

13cbeae

Signed-off-by: ErikQQY <2283984853@qq.com>

Fix k_interp

d79c04b

Signed-off-by: ErikQQY <2283984853@qq.com>

Merge branch 'master' into master

ebfd30e

ChrisRackauckas reviewed Jun 17, 2023

View reviewed changes

ErikQQY added 2 commits June 21, 2023 23:55

Fix MIRKTableau

6429528

Signed-off-by: ErikQQY <2283984853@qq.com>

Done mesh refinement

cb7ff37

Signed-off-by: ErikQQY <2283984853@qq.com>

Merge branch 'master' into master

3cbfecf

ErikQQY added 6 commits June 27, 2023 18:59

Refactor

1f59a70

Signed-off-by: ErikQQY <2283984853@qq.com>

Fix sum_stages

38d4eaa

Signed-off-by: ErikQQY <2283984853@qq.com>

Refactor MIRKGeneralCache

7fece70

Signed-off-by: ErikQQY <2283984853@qq.com>

Refactor k_interp

2409831

Signed-off-by: ErikQQY <2283984853@qq.com>

Fix BVPSystem in loop

89148bb

Signed-off-by: ErikQQY <2283984853@qq.com>

x in tableau should be s x s

711acd0

Signed-off-by: ErikQQY <2283984853@qq.com>

Fianlly working

7c2cc23

Signed-off-by: ErikQQY <2283984853@qq.com>

ErikQQY changed the title ~~[WIP] Add defect control adaptivity~~ Add defect control adaptivity Jul 10, 2023

avik-pal reviewed Jul 10, 2023

View reviewed changes

ErikQQY added 2 commits July 17, 2023 20:55

Refactor

dbebe2f

Signed-off-by: ErikQQY <2283984853@qq.com>

Refactor half_mesh

a13f1d4

Signed-off-by: ErikQQY <2283984853@qq.com>

ChrisRackauckas reviewed Jul 18, 2023

View reviewed changes

src/adaptivity.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Jul 18, 2023

View reviewed changes

src/adaptivity.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Jul 18, 2023

View reviewed changes

src/adaptivity.jl

Generate a new mesh.

"""

function redistribute(mesh_current::Vector,

n::Int64,

Copy link

Member

ChrisRackauckas Jul 18, 2023

Choose a reason for hiding this comment

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

Why require Int64?

ChrisRackauckas reviewed Jul 18, 2023

View reviewed changes

src/adaptivity.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Jul 18, 2023

View reviewed changes

src/adaptivity.jl Show resolved Hide resolved

ChrisRackauckas reviewed Jul 18, 2023

View reviewed changes

src/adaptivity.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Jul 18, 2023

View reviewed changes

src/adaptivity.jl Outdated Show resolved Hide resolved

ErikQQY added 2 commits July 21, 2023 20:45

Refactor

a02c6e0

Signed-off-by: ErikQQY <2283984853@qq.com>

Fix retcode of NonlinearSolve

14d5a34

Signed-off-by: ErikQQY <2283984853@qq.com>

ErikQQY and others added 2 commits July 23, 2023 11:30

Merge branch 'master' into master

1009899

All done

6a4860f

Signed-off-by: ErikQQY <2283984853@qq.com>

ChrisRackauckas merged commit e697350 into SciML:master Jul 24, 2023
4 of 6 checks passed

ErikQQY mentioned this pull request Jul 24, 2023

More TODOs #93

Open

4 tasks

sebastiantk mentioned this pull request Mar 18, 2024

MIRK4 adaptive time steps #172

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add defect control adaptivity #89

Add defect control adaptivity #89

ErikQQY commented May 29, 2023 •

edited

ChrisRackauckas Jun 17, 2023

ErikQQY Jun 20, 2023

ChrisRackauckas Jun 17, 2023

ErikQQY Jun 20, 2023

ChrisRackauckas Jun 17, 2023

ErikQQY Jun 20, 2023

ChrisRackauckas Jun 17, 2023

ErikQQY Jun 20, 2023

ChrisRackauckas commented Jun 24, 2023

ErikQQY commented Jun 24, 2023 •

edited

avik-pal Jul 10, 2023

ErikQQY Jul 10, 2023

ChrisRackauckas Jul 18, 2023

ChrisRackauckas Jul 18, 2023

ChrisRackauckas Jul 18, 2023

ChrisRackauckas Jul 18, 2023

ChrisRackauckas Jul 18, 2023

ChrisRackauckas Jul 18, 2023

ChrisRackauckas Jul 18, 2023

ChrisRackauckas Jul 18, 2023

ChrisRackauckas Jul 18, 2023

ChrisRackauckas Jul 18, 2023

ErikQQY Jul 18, 2023 •

edited

ChrisRackauckas Jul 18, 2023

ChrisRackauckas Jul 18, 2023

ChrisRackauckas commented Jul 22, 2023

ErikQQY commented Jul 24, 2023

ChrisRackauckas commented Jul 24, 2023

	mesh_new[2 * i] = (mesh_current[i + 1] + mesh_current[i]) / 2.0
	mesh_new[2 * i] = (mesh_current[i + 1] + mesh_current[i]) / 2

	weights_2, weights_2_prime = interp_weights(1.0 - tau_star, alg)
	weights_2, weights_2_prime = interp_weights(1 - tau_star, alg)

Add defect control adaptivity #89

Add defect control adaptivity #89

Conversation

ErikQQY commented May 29, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ChrisRackauckas commented Jun 24, 2023

ErikQQY commented Jun 24, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ErikQQY Jul 18, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ChrisRackauckas commented Jul 22, 2023

ErikQQY commented Jul 24, 2023

ChrisRackauckas commented Jul 24, 2023

ErikQQY commented May 29, 2023 •

edited

ErikQQY commented Jun 24, 2023 •

edited

ErikQQY Jul 18, 2023 •

edited