UpwindDifference Banded concretization doesn't fill memory

[hurricane007]

Hello all,
I am not sure if this is an error. However, when I use UpwindDifference the result is not what I wanted and changes every time. For example, if I want a 1st order upwind differencing with 2nd order precision, I got the results below:

julia> BandedMatrix(UpwindDifference(1, 2, 0.01, 9))
9×11 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
0.0 -8.20111e-313 1.09348e-312 -2.7337e-313 ⋅ … ⋅ ⋅ ⋅
0.0 0.0 -8.20111e-313 1.09348e-312 -2.7337e-313 ⋅ ⋅ ⋅
⋅ 0.0 0.0 -8.20111e-313 1.09348e-312 ⋅ ⋅ ⋅
⋅ ⋅ 0.0 0.0 -8.20111e-313 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.0 0.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 0.0 … -2.7337e-313 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 1.09348e-312 -2.7337e-313 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ -8.20111e-313 1.09348e-312 -2.7337e-313
⋅ ⋅ ⋅ ⋅ ⋅ -2.7337e-313 0.0 2.7337e-313

julia> BandedMatrix(UpwindDifference(1, 2, 0.01, 9))
9×11 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
0.0 -1.15603e-313 1.54137e-313 -3.85342e-314 ⋅ … ⋅ ⋅ ⋅
0.0 0.0 -0.0 0.0 -0.0 ⋅ ⋅ ⋅
⋅ 0.0 0.0 -1.15616e-313 1.54155e-313 ⋅ ⋅ ⋅
⋅ ⋅ 0.0 0.0 -1.15678e-313 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.0 0.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 0.0 … -3.85388e-314 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 3.41183e-313 -8.52958e-314 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ -0.0 0.0 -0.0
⋅ ⋅ ⋅ ⋅ ⋅ -3.85388e-314 0.0 3.85388e-314

[ChrisRackauckas]

It looks like it's making the right structure but not filling in the values for some reason. Those e-314 are coming from a similar call that is leaving undef values. I think there's a chance a generic fallback is being hit.

I've been struggling with getting the KdV test working, and it seems to be related to this. I'll take a look.

It seems like this is caused by uninitialized data leaking into the coefficients attribute of UpwindDifference operators. Culprits here and (repeated) here, since compute_coeffs! reads uninitialized values from coefficients.

The same can happen with CenteredDifference if the coeff_func uses the current coefficients (see this line).

The tests in upwind_operators_interface.jl should catch this, but it seems to be an unlikely bug to surface for small grid sizes. I only regularly see it with N > 250.

Unless I'm mistaken, some code still assumes coeff_func is always a function (like here). I also don't understand why the first parameter is named du even though the value is treated as coefficients.

I'm tempted to simplify the interface to just having the coeff_func::Function case, and rely on the caller to create constant / spatially dependent functions. What do you think @ChrisRackauckas?

Also, if the plan is to eventually redo this with ModelingToolkit, that would be helpful to know.

[ChrisRackauckas]

We should probably initialize it to 1 at least: initializing to undef is indeed dangerous. I don't think there's an issue with allowing coefficients as vectors and such, because 5A is a very natural thing to have work, but since we do that, we should interpret A == 1*A.

[hurricane007]

Hi Chris and ghost,

there seems another issue when a function is used as the coeff_fun. I tried the example in UpwindDifference:

drift = [1., 1., -1.]
UpwindDifference(1, 1, 1., 3, i -> drift[i])

and the output is :

BoundsError: attempt to access 3-element Array{Float64,1} at index [[5.65079456e-315, 5.65079472e-315, 5.650794877e-315]]

Stacktrace:
[1] throw_boundserror(::Array{Float64,1}, ::Tuple{Array{Float64,1}}) at .\abstractarray.jl:537
[2] checkbounds at .\abstractarray.jl:502 [inlined]
[3] _getindex at .\multidimensional.jl:726 [inlined]
[4] getindex(::Array{Float64,1}, ::Array{Float64,1}) at .\abstractarray.jl:980
[5] (::var"#7#8")(::Array{Float64,1}) at .\In[10]:1
[6] compute_coeffs!(::var"#7#8", ::Array{Float64,1}) at C:\ProgramData.juia\packages\DiffEqOperators\QiB3a\src\derivative_operators\coefficient_functions.jl:18
[7] UpwindDifference{1}(::Int64, ::Int64, ::Float64, ::Int64, ::Function) at C:\ProgramData.juia\packages\DiffEqOperators\QiB3a\src\derivative_operators\derivative_operator.jl:174
[8] UpwindDifference(::Int64, ::Vararg{Any,N} where N) at C:\ProgramData.juia\packages\DiffEqOperators\QiB3a\src\derivative_operators\derivative_operator.jl:236
[9] top-level scope at In[10]:1

so the problem is in coefficient_function.jl line 18:
current_coeffs[:] = coeff_func(current_coeffs)
where current_coeffs was initialized as
Vector{T}(undef,len)
in the example above, the function coeff_func accepts the index number.
I want to modify this function. There is, however, another possibility that the function returns the coefficient according to its coordinate, e.g., coeff_fun is 'x -> x^2 + 1. '
I am not really into Julia. So I have a question: is the matrix generated by UpWindDifference function variable? Because in, e.g., CFD, the local velocity may change with time. Thus the direction of UpWind changes.
Sorry I also don't know how the UpWind scheme is implemented in real CFD code. I am not able to modify this function. If I have some new idea I will come back.

BTW, i->drift[i] is used as the function for CenteredDifference, there is another issue as well. With the same definition of drift,

BandedMatrix(CenteredDifference(1, 2, 0.1, 5, drift))
gives

3×5 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
-9.10811e-315 0.0 9.10811e-315 0.0 ⋅
0.0 -9.10811e-315 0.0 9.10811e-315 0.0
0.0 0.0 -9.10812e-315 0.0 9.10812e-315
which also looks like a undefined matrix. But seems the centered difference has nothing to do with the symbol of the coefficient, yes?

Regards,

[ChrisRackauckas]

o I have a question: is the matrix generated by UpWindDifference function variable? Because in, e.g., CFD, the local velocity may change with time. Thus the direction of UpWind changes.

It's not a matrix so that it can switch directions easily. Just modify the coefficients and it'll change its action.

[hurricane007]

thanks for fixing the issue. However, I found an issue concerning UpwindDifference when the coeff_fun is not given. For example:
julia> drift=[1.0, 1., -1., 1., -1.]
julia> BandedMatrix(UpwindDifference(1,1,0.1,5,drift))
gives correct result.

5×7 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
0.0 -10.0 10.0 ⋅ ⋅ ⋅ ⋅
⋅ 0.0 -10.0 10.0 ⋅ ⋅ ⋅
⋅ ⋅ 10.0 -10.0 0.0 ⋅ ⋅
⋅ ⋅ ⋅ 0.0 -10.0 10.0 ⋅
⋅ ⋅ ⋅ ⋅ 10.0 -10.0 0.0

however,
julia> BandedMatrix(UpwindDifference(1,1,0.1,5))

gives wrong result

5×7 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
0.0 -1.93627e-314 1.93627e-314 ⋅ ⋅ ⋅ ⋅
⋅ 0.0 -1.93627e-314 1.93627e-314 ⋅ ⋅ ⋅
⋅ ⋅ 0.0 -1.93627e-314 1.93627e-314 ⋅ ⋅
⋅ ⋅ ⋅ 0.0 -1.93627e-314 1.93627e-314 ⋅
⋅ ⋅ ⋅ ⋅ 0.0 -1.93627e-314 1.93627e-314

maybe UpwindDifference without coeff_func should not be allowed?
and the problem concerning CenteredDifference remains,
julia> Array(CenteredDifference(1,2,0.1, 5,drift))

5×7 Array{Float64,2}:
-4.09934e-315 0.0 4.09934e-315 0.0 0.0 0.0 0.0
0.0 -4.03688e-315 0.0 4.03688e-315 0.0 0.0 0.0
0.0 0.0 -4.16402e-315 0.0 4.16402e-315 0.0 0.0
0.0 0.0 0.0 -3.99142e-315 0.0 3.99142e-315 0.0
0.0 0.0 0.0 0.0 -3.99142e-315 0.0 3.99142e-315