Rebased #168

[ChrisRackauckas]

No description provided.

[ChrisRackauckas]

@grahamas @xtalax when you rebase you have to force push, not do a merge afterwards. So I cherrypicked one commit back and it's all good.

[ChrisRackauckas]

how did these break?

[grahamas]

edit: disregard this comment

one of my "rebase" commits changes A.L.len to size(A.L, 2) in the definition of size(::GhostDerivativeOperator) which all the broken tests rely on. I don't see any reason whatsoever why I would have done that -- both master and xtalax have A.L.len and I'm aware that those aren't the same thing. Best guess is I was being stupid. So I can fix that, but where should I commit?

[ChrisRackauckas]

PR to this branch.

[grahamas]

Unfortunately I can't fix this. In doing this rebase, I was purely fixing things that didn't require me to make changes to what I saw as @xtalax's intentions/the code's semantics.

My mistake was today, not in the rebase. Today, I was mistakenly looking at xtalax:master rather than xtalax:#168 to figure out what the outcome is supposed to be.

In fact, xtalax:#168 does redefine size(A::GhostDerivativeOperator) = (size(A.L, 2),size(A.L, 2)). But the constructor for CenteredDifference(_, _, 10) will very definitely have second dimension of size 10+2. I don't know which definition of size (or maybe the constructor?) is wrong. So two different parts of the code have two clear but conflicting effects. I can start trying to fix the code's semantics, but that'll be a longer process.

edit: concisely, branch #168 breaks this test by redefining size(::GhostDerivativeOperator) and the rebase faithfully also breaks the test.

[ChrisRackauckas]

I would revert that change I think?

[grahamas]

Reverted and passing. @ChrisRackauckas how often do you want me to PR from my personal branch into this branch?

[ChrisRackauckas]

@xtalax note that this isn't needed anymore now that we are using @safetestsets

[xtalax]

noted, i'll change that around on my next commit

[ChrisRackauckas]

Not sure if this is perfect, but perfect might be the enemy of good here. Let's make sure this is passing tests and get it merged. Then we can continue to clean up the other little pieces. I think the next thing I want to do is to try and start getting a real documentation together for this library to really take stock of where we're at.

[ChrisRackauckas]

It looks like there's a lot of failing tests? Are there commits that were not on that branch?

[ChrisRackauckas]

Doesn't look like this can merge because there's a ton of broken tests.

[grahamas]

Something's weird -- I'm not sure this is from the right commit. Line 23 of test/runtests.jl is causing an error that's fixed post-rebase (cf. https://github.com/xtalax/DiffEqOperators.jl/pull/26/files#diff-fce720c43af3c52c862fd7451c7374b8 which shows how I'd expect this file to look)

[grahamas]

That's my guess, too, but I don't know why @xtalax changed it in the first place, or if other new code relies on the new definition. If this PR is still unresolved at the end of June (after a thesis committee meeting) I'd be happy to jump in completely. Until then I wouldn't be comfortable asserting a change like this. Still happy to resolve any rebase-related problems in the meantime, since those are more mechanical and local. Though depending on how long it takes to resolve the failing tests, I might just need to run another rebase in a month....

…

On Sun, Jun 7, 2020, 9:56 AM Christopher Rackauckas < ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In test/ghost_derivative_operators_interface.jl <#246 (comment)> : > @@ -7,7 +7,7 @@ using DiffEqOperators, Test A = L * Q @test ndims(A) == 2 - @test size(A) == (10, 10) - @test size(A, 1) == 10 - @test length(A) == 100 + @test_broken size(A) == (10, 10) I would revert that change I think? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#246 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AANS24FXNTTS3IGWBTQTIQ3RVOTCBANCNFSM4NVO4GSQ> .

[ChrisRackauckas]

No worries. I think these test failures are what stalled it in the first place, so we can work through them whenever there's time.

[xtalax]

@ChrisRackauckas this error and those just below are what stalled progress before - I really have no idea why this isn't working

[ChrisRackauckas]

what's the error?

[xtalax]

the test marked broken should definitely pass if the above three are passing, and yet it doesn't

[grahamas]

It feels like it has to be a bug in the one of the \ implementations (or, more likely, expected imprecision). Running everything in the REPL:

EDIT: added fourth_deriv_approx_stencil definition; other things that I accidentally got from including test file

julia> using DiffEqOperators

julia> function fourth_deriv_approx_stencil(N)
    A = zeros(N,N+2)
    A[1,1:8] = [3.5 -56/3 42.5 -54.0 251/6 -20.0 5.5 -2/3]
    A[2,1:8] = [2/3 -11/6 0.0 31/6 -22/3 4.5 -4/3 1/6]

    A[N-1,N-5:end] = reverse([2/3 -11/6 0.0 31/6 -22/3 4.5 -4/3 1/6], dims=2)
    A[N,N-5:end] = reverse([3.5 -56/3 42.5 -54.0 251/6 -20.0 5.5 -2/3], dims=2)

    for i in 3:N-2
        A[i,i-2:i+4] = [-1/6 2.0 -13/2 28/3 -13/2 2.0 -1/6]
    end
    return A
end

julia> dx = 0.01; x = 0.01:dx:0.2; N=length(x); u=sin.(x); L = CenteredDifference(4, 4, dx, N); Q = RobinBC((1.0, 0.0, 0.0), (1.0, 0.0, sin(0.2+dx)), dx); A = L*Q;                         
                                                                                              
julia> analytic_L = fourth_deriv_approx_stencil(N) ./ dx^4; analytic_QL = [transpose(zeros(N)); Diagonal(ones(N)); transpose(zeros(N))]; analytic_AL = analytic_L * analytic_QL; analytic_Qb = [zeros(N+1); sin(0.2+dx)]; analytic_Ab = analytic_L*analytic_Qb; 

julia> using SparseArrays; As_l, As_b = sparse(A, length(u))

julia> analytic_AL \ Vector(u .- analytic_Ab) ≈ As_l \ Vector(u .- As_b)
false

julia> analytic_AL \ Vector(u .- analytic_Ab) ≈ Array(As_l) \ Vector(u .- As_b)
true

julia> methods(\, typeof.([analytic_AL, Vector(u .- analytic_Ab)]))
# 1 method for generic function "\":                                                                                         
[1] \(A::AbstractArray{T,2} where T, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T)
 in LinearAlgebra at /home/graham/.local/opt/julia-1.5.0/share/julia/stdlib/v1.5/LinearAlgebra
/src/generic.jl:1102

julia> methods(\, typeof.([As_l, Vector(u .- As_b)]))
# 1 method for generic function "\":
[1] \(A::SparseArrays.AbstractSparseMatrixCSC, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) in SparseArrays at /home/graham/.local/opt/julia-1.5.0/share/julia/stdlib/v1.5/SparseArrays/src/linalg.jl:1460  

[grahamas]

I filed an issue with an MWE: JuliaLang/julia#37004. I'll mark this as broken for now and move on.

[grahamas]

They pointed out that the difference is within the expected range given the condition number of analytic_AL. I've revised the test to reflect that, and now it and a nearby similar test are passing. I imagine this is a problem in some other tests as well. Would we rather use a better conditioned matrix? I wouldn't know how to choose.

[ChrisRackauckas]

Just change the tolerance to be fit for the condition number. Any stencil operation is going to be pretty ill-conditioned (this is where the stiffness comes from)

[grahamas]

Status update: Only two problem areas remaining.

1. upwind_operators_interface.jl (test file) produces errors at random, varying between runs. This obviously drives me insane, especially since I'm not even sure how the ghost derivative stuff should be affecting upwind operators, so I'm leaving it for last.

2. 3D_laplacian.jl (test file) reveals that there's no mul! method for multidimensional operators -- all existing mul! methods assert that the padding is in only one dimension. A cursory dive back into the pre-rebase commit history suggests I didn't accidentally delete anything. That being the case, I'll have to rework the fallback in derivative_operator_functions.jl to do a check for higher dimensional padding.

[grahamas]

The fallback mul! now works with multidimensionally padded inputs by introducing unpadded_size helper function. TODO:

1. Update specialized mul!s that call NNlib.conv! [done!]
2. R̶e̶d̶e̶f̶i̶n̶e̶ ̶̶A̶*̶Q̶̶ ̶f̶o̶r̶ ̶̶A̶:̶:̶D̶i̶f̶f̶E̶q̶O̶p̶e̶r̶a̶t̶o̶r̶C̶o̶m̶b̶i̶n̶a̶t̶i̶o̶n̶̶ ̶(̶j̶u̶s̶t̶ ̶̶A̶<̶:̶A̶b̶s̶t̶r̶a̶c̶t̶D̶i̶f̶f̶E̶q̶O̶p̶e̶r̶a̶t̶o̶r̶̶?̶)̶ ̶a̶n̶d̶ ̶̶Q̶:̶:̶C̶o̶m̶p̶o̶s̶e̶d̶M̶u̶l̶t̶i̶D̶i̶m̶B̶C̶̶ ̶(̶j̶u̶s̶t̶ ̶̶Q̶<̶:̶M̶u̶l̶t̶i̶D̶i̶m̶B̶C̶̶?̶)̶ ̶a̶s̶ ̶a̶ ̶c̶o̶m̶p̶o̶s̶i̶t̶i̶o̶n̶ ̶s̶o̶ ̶t̶h̶a̶t̶ ̶̶A̶*̶Q̶*̶u̶̶ ̶j̶u̶s̶t̶ ̶w̶o̶r̶k̶s̶ ̶(̶a̶l̶t̶e̶r̶n̶a̶t̶i̶v̶e̶ ̶i̶s̶ ̶r̶e̶p̶a̶d̶d̶i̶n̶g̶ ̶̶u̶̶ ̶f̶o̶r̶ ̶e̶v̶e̶r̶y̶ ̶d̶i̶m̶e̶n̶s̶i̶o̶n̶.̶.̶.̶ ̶s̶e̶e̶m̶s̶ ̶u̶n̶w̶i̶s̶e̶)̶ hah I almost tried to reinvent GhostDerivativeOperators -- caught myself just in time.

[ChrisRackauckas]

Awesome. Should this be passing then?

[ChrisRackauckas]

Looks like there's a merge conflict.

[grahamas]

@ChrisRackauckas Nope! I don't have access to this branch; I was working on yet another fork. I've opened #279 which is rebased to master as of 20 minutes ago and is passing all tests on my machine.

[ChrisRackauckas]

alright, let's continue there.

[ChrisRackauckas]

Let's zero this so errors with undef memory don't pop up.