add J, W cache and linsolve interface for rosenbrock oop methods and add J cache for firk oop #873

huanglangwen · 2019-08-14T09:09:53Z

devmotion · 2019-08-14T09:25:25Z

I think we should not introduce an additional implementation for generating J (and W) but rather use/modify the existing
https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/cd45039e49eefdc5d60aa72b54f91e84851eb747/src/integrators/integrator_utils.jl#L386
and
https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/cd45039e49eefdc5d60aa72b54f91e84851eb747/src/integrators/integrator_utils.jl#L412

codecov · 2019-08-14T10:04:40Z

Codecov Report

Merging #873 into master will increase coverage by 0.18%.
The diff coverage is 85.05%.

@@            Coverage Diff             @@
##           master     #873      +/-   ##
==========================================
+ Coverage   80.16%   80.35%   +0.18%     
==========================================
  Files          92       92              
  Lines       30269    30244      -25     
==========================================
+ Hits        24265    24302      +37     
+ Misses       6004     5942      -62

Impacted Files	Coverage Δ
src/caches/rosenbrock_caches.jl	`99.58% <100%> (+14.77%)`	⬆️
src/tableaus/rosenbrock_tableaus.jl	`99.19% <100%> (+0.02%)`	⬆️
src/generic_rosenbrock.jl	`98.3% <100%> (ø)`	⬆️
src/caches/firk_caches.jl	`93.75% <50%> (-2.81%)`	⬇️
src/integrators/integrator_utils.jl	`86.72% <65.62%> (+1.55%)`	⬆️
src/caches/adams_bashforth_moulton_caches.jl	`80.7% <0%> (+0.44%)`	⬆️
src/caches/kencarp_kvaerno_caches.jl	`97.02% <0%> (+1.98%)`	⬆️
src/caches/sdirk_caches.jl	`96.52% <0%> (+3.47%)`	⬆️
src/caches/bdf_caches.jl	`89.54% <0%> (+3.92%)`	⬆️
... and 1 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update cd45039...dee2e87. Read the comment docs.

huanglangwen · 2019-08-14T13:51:25Z

It seems the oopgenW is broken when u is a StaticArray. Not sure how to fix it.
https://travis-ci.org/JuliaDiffEq/OrdinaryDiffEq.jl/jobs/571800687

devmotion · 2019-08-14T14:07:01Z

I guess it's the first time that someone actually tests the whole logic of oop_generate_W. I think
https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/cd45039e49eefdc5d60aa72b54f91e84851eb747/src/integrators/integrator_utils.jl#L426-L432
is incorrect, since column vectors with only one entry end up in the second branch. I assume it should be changed to

J = if ndims(u) == 2 # `u` is already a matrix
        u
      else # `u` is a vector
        hcat(u, u)
      end

Can you check if that solves the issue?

To be honest, it still feels very crude to me, the units are incorrect and will fail if the dimension of u is greater than 2.

ChrisRackauckas · 2019-08-14T14:09:14Z

Aren't we already vecing u when applying the linear operator? If not, we should, and in that case false .* rate_prototype .* rate_prototype' should work, ever for u a Number.

devmotion · 2019-08-14T14:19:05Z

Aren't we already vecing u when applying the linear operator? If not, we should, and in that case false .* rate_prototype .* rate_prototype' should work, ever for u a Number.

I don't see how this is related to the test failure? But I agree that probably we should just replace the snippet above and https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/cd45039e49eefdc5d60aa72b54f91e84851eb747/src/integrators/integrator_utils.jl#L438 with

J = false .* _vec(u) .* _vec(u)'

That would be much cleaner.

Then we could just branch depending on the type of u for the computation of W (if we do not want to compute one LU factorization, see #672).

ChrisRackauckas · 2019-08-14T14:21:18Z

The real difficulty is that we don't have an OOP linsolve interface, which we need to get in order to really complete this I think?

devmotion · 2019-08-14T14:22:00Z

And the fact that lu(J) will fail since J is singular, as I just noticed...

huanglangwen · 2019-08-14T14:25:54Z

And the fact that lu(J) will fail since J is singular, as I just noticed...

Your patch works anyway. Shall I push that?

devmotion · 2019-08-14T14:39:47Z

And the fact that lu(J) will fail since J is singular, as I just noticed...

Actually, it seems lu(@SMatrix zeros(2,2)) works (and only lu(zeros(2,2)) fails), so I think you could even write

J = false .* _vec(u) .* _vec(u)'
W = if u isa StaticArray
  lu(J)
  elseif u isa Number
    u
  else
    LU{LinearAlgebra.lutype(uEltypeNoUnits)}(Matrix{uEltypeNoUnits}(undef, 0, 0),
                                                 Vector{LinearAlgebra.BlasInt}(undef, 0),
                                                 zero(LinearAlgebra.BlasInt))
  end

instead of https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/1330fa7a01626a4c334d171a619802475be8d036/src/integrators/integrator_utils.jl#L432-L448

huanglangwen · 2019-08-14T14:46:30Z

The real difficulty is that we don't have an OOP linsolve interface, which we need to get in order to really complete this I think?

What's the motivation to move from k₃ = _reshape(W\-_vec(linsolve_tmp), axes(uprev))?

YingboMa · 2019-08-14T16:53:09Z

W\b is wasteful to do at every stage. A better strategy is to factorize it first and reuse the factorized object. Hence we need to add an OOP linsolve interface.

YingboMa · 2019-08-14T16:55:04Z

~~W\b is wasteful to do at every stage. A better strategy is to factorize it first and reuse the factorized object.~~

Never mind. We want to use factorization methods other than lu, so we need it.

src/integrators/integrator_utils.jl

Co-Authored-By: Yingbo Ma <[email protected]>

devmotion · 2019-08-15T06:01:02Z

src/caches/rosenbrock_caches.jl

  tmp = zero(rate_prototype)
  atmp = similar(u, uEltypeNoUnits)
-  tab = Rosenbrock23ConstantCache(real(uBottomEltypeNoUnits),identity,identity)
+  tab = Rosenbrock23ConstantCache(real(uBottomEltypeNoUnits),identity,identity,nothing,nothing,nothing)


I always thought that we shouldn't misuse the oop cache as tableau but rather use separate tableau structs - and even more so with the new fields for J and W.

yeah, it's really just a short hand that always worked... until it didn't.

huanglangwen · 2019-08-15T12:06:23Z

Is this pr good for merge? I want to start a new pr for unifying oop and iip in derivative_utils.jl .

devmotion · 2019-08-15T12:08:05Z

src/caches/rosenbrock_caches.jl

  tmp = zero(rate_prototype)
  atmp = similar(u, uEltypeNoUnits)
-  tab = Rosenbrock23ConstantCache(real(uBottomEltypeNoUnits),identity,identity,nothing,nothing,nothing)
+  tab = Ros23ConstantCache(real(uBottomEltypeNoUnits))


I think I would prefer the name Rosenbrock23Tableau or something similar, Ros23ConstantCache seems a bit strange to me.

I agree, but then we need to correct all the tab struct in \tableaus as well as places calling the tabs in \caches

Sure, there might be some other structs which probably should be renamed as well but that can be done in a different PR, I guess. For now, IMO it's sufficient to name the new struct in a reasonable way 😃

devmotion · 2019-08-15T12:10:12Z

src/caches/rosenbrock_caches.jl

  tmp = zero(rate_prototype)
  atmp = similar(u, uEltypeNoUnits)
-  tab = Rosenbrock32ConstantCache(real(uBottomEltypeNoUnits),identity,identity,nothing,nothing,nothing)
+  tab = Ros32ConstantCache(real(uBottomEltypeNoUnits))


Similar to my comment above, I think we should use something a bit more meaningful than Ros32ConstantCache.

ChrisRackauckas · 2019-08-16T04:42:23Z

Let's continue the renames as another PR. Also note that in the sphere of influence for this should be @saurabhkgp21's newest implicit extrapolation methods. For those methods, we need a safe-hatch so that way Jacobian reuse is just generally turned off. In other cases, having an easy switch on it is nice anyways, so it's something to keep in mind.

add J and W cache for rosenbrock and firk oop

c4b6c93

refactor init J and W

c2d9f45

fix test error

1330fa7

fix test error

9a46866

YingboMa reviewed Aug 14, 2019

View reviewed changes

src/integrators/integrator_utils.jl Outdated Show resolved Hide resolved

huanglangwen and others added 2 commits August 15, 2019 11:07

Update src/integrators/integrator_utils.jl

1653cb9

Co-Authored-By: Yingbo Ma <[email protected]>

add linsolve interface to oop rosenbrock methods

a59f756

huanglangwen changed the title ~~add J and W cache for rosenbrock and firk oop~~ add J, W cache and linsolve interface for rosenbrock oop methods and add J cache for firk oop Aug 15, 2019

devmotion reviewed Aug 15, 2019

View reviewed changes

split tab and oop cache for ros23/32

108d091

devmotion reviewed Aug 15, 2019

View reviewed changes

#TODO: Rename other tableau structs

dee2e87

ChrisRackauckas merged commit 893bc16 into SciML:master Aug 16, 2019

huanglangwen mentioned this pull request Aug 16, 2019

[WIP] Simplify derivative_utils.jl #876

Merged

huanglangwen mentioned this pull request Aug 16, 2019

Rename tableau struct for rosenbrock methods #877

Merged

Uh oh!

add J, W cache and linsolve interface for rosenbrock oop methods and add J cache for firk oop #873

add J, W cache and linsolve interface for rosenbrock oop methods and add J cache for firk oop #873

Uh oh!

Conversation

huanglangwen commented Aug 14, 2019

Uh oh!

devmotion commented Aug 14, 2019

Uh oh!

codecov bot commented Aug 14, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

huanglangwen commented Aug 14, 2019

Uh oh!

devmotion commented Aug 14, 2019

Uh oh!

ChrisRackauckas commented Aug 14, 2019

Uh oh!

devmotion commented Aug 14, 2019

Uh oh!

ChrisRackauckas commented Aug 14, 2019

Uh oh!

devmotion commented Aug 14, 2019

Uh oh!

huanglangwen commented Aug 14, 2019

Uh oh!

devmotion commented Aug 14, 2019

Uh oh!

huanglangwen commented Aug 14, 2019

Uh oh!

YingboMa commented Aug 14, 2019

Uh oh!

YingboMa commented Aug 14, 2019

Uh oh!

Uh oh!

devmotion Aug 15, 2019

Choose a reason for hiding this comment

Uh oh!

ChrisRackauckas Aug 15, 2019

Choose a reason for hiding this comment

Uh oh!

huanglangwen commented Aug 15, 2019

Uh oh!

devmotion Aug 15, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

huanglangwen Aug 15, 2019

Choose a reason for hiding this comment

Uh oh!

devmotion Aug 15, 2019

Choose a reason for hiding this comment

Uh oh!

devmotion Aug 15, 2019

Choose a reason for hiding this comment

Uh oh!

ChrisRackauckas commented Aug 16, 2019

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

codecov bot commented Aug 14, 2019 •

edited

Loading

devmotion Aug 15, 2019 •

edited

Loading