NaN when using space grid with x_0 = 0

[hurricane007]

Hello, I am not sure if it's an issue. But when I using a space grid with x0 = 0, the first two elements of UpwindDifference(1,2, Δx, length(x), -usg)bcu0 are NaN, and the first elements of CenteredDifference(1,2,Δx, length(x))bcu0 is NaN. Is it an error? or the spac discretization must start from x0 > 0?

Test code:

using DiffEqOperators
Δx = 0.1
x0 = 0.0
x = x0 :Δx:1.
u0 = exp.(1 .-x) .-1 
Deg = 1E-4 # m^2/s
usg = 0.5
bc = RobinBC([1., -Deg/usg, 1. ], [0. , 1 , 0], x,1.)
UpwindDifference(1,2, Δx, length(x), -usg)*bc*u0

11-element Array{Float64,1}:
 NaN
 NaN
   1.1087695767283758
   1.0032562010037278
   0.9077837505447812
   0.8213967049779365
   0.7432304737154816
   0.6725027428423593
   0.608505645455582
   0.5505986770943327
   0.4982022853560525

and
CenteredDifference(1,2,Δx, length(x))*bc*u0

11-element Array{Float64,1}:
 NaN
  -2.463704499832886
  -2.2292520184323656
  -2.0171106405097943
  -1.8251571838517426
  -1.651470513746193
  -1.494312315620625
  -1.3521096974055025
  -1.2234394475017774
  -1.1070137908008493
  -0.5258545903782386

[hurricane007]

In addition, this tow also give different result, does it mean the "ghost node" is a real node at x = 0?

UpwindDifference(1,2, Δx, length(range(0.01, step=Δx, length=11)), -usg)*
RobinBC([1., -Deg/usg, 1.2], [0. , 1 , 0], (range(0.01, step=Δx, length=11)),1.)*
(exp.(1 .-(range(0.01, step=Δx, length=11))) .-1)

UpwindDifference(1,2, Δx, length(range(0.1, step=Δx, length=11)), -usg)*
RobinBC([1., -Deg/usg, 1.2], [0. , 1 , 0], (range(0.1, step=Δx, length=11)),1.)*
(exp.(1 .-(range(0.1, step=Δx, length=11))) .-1)

[hurricane007]

Beside, if I want to solve a equation:
du/dt = -usg*du/dz + Deg d^2u/dz^2 + r(u)
is the defination below correct?

D1 = UpwindDifference(1,2, Δx, length(range(0.01, step=Δx, length=11)), -usg)
D2 = CenteredDifference(2,2,Δx, length(x))
bc = RobinBC([1., -Deg/usg, 1. ], [0. , 1 , 0], x,1.)
pdeFun(u,p,t) = D1*bc*u + Deg*bc*u + r(u)

[ChrisRackauckas]

That looks right. This looks like a bug, and it's not immediately clear where it comes from but it's worth looking into. I know we have a bunch of examples using x_0=0 so that's not it, it's something a little deeper.

[hurricane007]

Hello Chris, an update on this issue. I used x = x0 :Δx:1. in this example and get an error. But when I use x = collect(x0 :Δx:1.), the error disappeared. So you may know where the problem is.

[ChrisRackauckas]

Ahh thanks. I don't know where the bug is from that but it helps a ton. I don't know if I'll get to this "soon" so do use the workaround, but I hope to make sure this gets tackled as part of a summer project soon.

[hurricane007]

Ahh thanks. I don't know where the bug is from that but it helps a ton. I don't know if I'll get to this "soon" so do use the workaround, but I hope to make sure this gets tackled as part of a summer project soon.

Hi Chris, I am back again for this issue lol.
I didn't have time to continue the work with Julia until today.
And yes, there is still a problem when I start from 0, e.g. , in the example heat equation, https://github.com/SciML/DiffEqOperators.jl/blob/master/examples/heat_equation.jl , when I change boundary condition to

knots = collect(range(0.0, step=h, length=nknots))
bc = RobinBC([0.1, -0.001, 0.3], [0.0, 1.0, 0.0],knots)

the solution is totally different than using

knots = collect(range(h, step=h, length=nknots)); 
bc = RobinBC([0.1, -0.001, 0.3], [0.0, 1.0, 0.0],knots)

Seems the result indicates that the space discretization cannot start from 0.0. Actually, the generated boundary condition is different when starting spacing discretization at 0.0 or other.

knots = collect(range(0.0, step=h, length=nknots))
bc = RobinBC([0.1, -0.001, 0.3], [0.0, 1.0, 0.0],knots)
RobinBC{Float64,Array{Float64,1}}([0.0], 0.0, [1.0102040816326532], 0.0)

but

knots = collect(range(h, step=h, length=nknots));
bc = RobinBC([0.1, -0.001, 0.3], [0.0, 1.0, 0.0],knots)
RobinBC{Float64,Array{Float64,1}}([0.2512437810945274], 1.4925373134328357, [1.0101010101010102], 0.0)

Looks I need to write code with boundary condition and discretization...

[ChrisRackauckas]

That is definitely odd... thanks for reporting.

[mjsheikh]

Hello, I am not sure if it's an issue. But when I using a space grid with x0 = 0, the first two elements of UpwindDifference(1,2, Δx, length(x), -usg)_bc_u0 are NaN, and the first elements of CenteredDifference(1,2,Δx, length(x))_bc_u0 is NaN. Is it an error? or the spac discretization must start from x0 > 0?

Test code:

using DiffEqOperators
Δx = 0.1
x0 = 0.0
x = x0 :Δx:1.
u0 = exp.(1 .-x) .-1 
Deg = 1E-4 # m^2/s
usg = 0.5
bc = RobinBC([1., -Deg/usg, 1. ], [0. , 1 , 0], x,1.)
UpwindDifference(1,2, Δx, length(x), -usg)*bc*u0

11-element Array{Float64,1}:
 NaN
 NaN
   1.1087695767283758
   1.0032562010037278
   0.9077837505447812
   0.8213967049779365
   0.7432304737154816
   0.6725027428423593
   0.608505645455582
   0.5505986770943327
   0.4982022853560525

and
CenteredDifference(1,2,Δx, length(x))*bc*u0

11-element Array{Float64,1}:
 NaN
  -2.463704499832886
  -2.2292520184323656
  -2.0171106405097943
  -1.8251571838517426
  -1.651470513746193
  -1.494312315620625
  -1.3521096974055025
  -1.2234394475017774
  -1.1070137908008493
  -0.5258545903782386

@hurricane007 Please note the following :

• RobinBC takes the step size vector as argument and seems you are passing the domain grid points instead. Please try passing Δx = 0.1*ones(10) in place of x.
• u0 comprises of interior points, you have also included the values at ghost nodes 0 and 1 as well.
• In UpwindDifference(1,2, Δx, length(x), -usg)bcu0 , it should be length(Δx) - 1 instead.

I hope this resolves your issues.

[hurricane007]

Hello, I am not sure if it's an issue. But when I using a space grid with x0 = 0, the first two elements of UpwindDifference(1,2, Δx, length(x), -usg)_bc_u0 are NaN, and the first elements of CenteredDifference(1,2,Δx, length(x))_bc_u0 is NaN. Is it an error? or the spac discretization must start from x0 > 0?
Test code:

using DiffEqOperators
Δx = 0.1
x0 = 0.0
x = x0 :Δx:1.
u0 = exp.(1 .-x) .-1 
Deg = 1E-4 # m^2/s
usg = 0.5
bc = RobinBC([1., -Deg/usg, 1. ], [0. , 1 , 0], x,1.)
UpwindDifference(1,2, Δx, length(x), -usg)*bc*u0

11-element Array{Float64,1}:
 NaN
 NaN
   1.1087695767283758
   1.0032562010037278
   0.9077837505447812
   0.8213967049779365
   0.7432304737154816
   0.6725027428423593
   0.608505645455582
   0.5505986770943327
   0.4982022853560525

and
CenteredDifference(1,2,Δx, length(x))*bc*u0

11-element Array{Float64,1}:
 NaN
  -2.463704499832886
  -2.2292520184323656
  -2.0171106405097943
  -1.8251571838517426
  -1.651470513746193
  -1.494312315620625
  -1.3521096974055025
  -1.2234394475017774
  -1.1070137908008493
  -0.5258545903782386

@hurricane007 Please note the following :

• RobinBC takes the step size vector as argument and seems you are passing the domain grid points instead. Please try passing Δx = 0.1*ones(10) in place of x.
• u0 comprises of interior points, you have also included the values at ghost nodes 0 and 1 as well.
• In UpwindDifference(1,2, Δx, length(x), -usg)_bc_u0 , it should be length(Δx) - 1 instead.

I hope this resolves your issues.

thanks for the reply! This really helped!