"Upwind Operators"

[sciml]

"Upwind Operators"

[ChrisRackauckas]

we should supply the required information as type parameters as the direction of wave is unlikely to change frequently

No, the direction will change frequently. This is pretty much what makes the problem difficult. In a nonlinear problem, this is dependent on whatever multiplies the u_x term. In KdV, this is easy since 6u*u_x will always be the same direction on the chosen interval. But do something like (u-3)*u_x and now the most important term for the evolution does change for maybe a few points every timestep.

So direction needs to be mutable. I'm thinking something like an array of Bools (or maybe a BitArray), interpreting true as forward and false as backward. Every time step, a user would need to update_directions!(A,coefficient), where you check the sign.(x) on the terms of the terms in the coefficient.

Now that doesn't cover every case. Think about this one:

https://en.wikipedia.org/wiki/Upwind_differencing_scheme_for_convection

In this case, the Peclet number (realted to that flow rate, but divide out diffusion) determines if upwind differencing is stable. So in this case, we have three regimes:

1. Direction is positive and large enough, upwind.
2. Direction is negative and large enough, downwind.
3. Direction is not large enough, central.

So maybe it's best to take (-1,0,1) from the user as possible values for update_directions!(A,coefficient), so they can denote the cutoffs themselves.

where +0.7 tells us the fraction of points in the stencil for a point which lie in front of it. So automatically 30% of points lie behind it.

This won't work. The key is not the percentage of points which need up/downwind, its which points need upwind/downwind. It's very common to have a few domains [x_i,x_i + dx_i] (of changing location and size)` which are downwind, in the middle of an operator where most everything else is upwind (eddies in fluid flows). It's important to know which domains are upwind and which domains are downwind.

[ChrisRackauckas]

This is a good paper which describes pretty much what we want in the end.

http://www.sciencedirect.com/science/article/pii/S0377042705000634#bib5

[shivin9]

@ChrisRackauckas @dextorious Any thoughts on the API for integrating this?

[shivin9]

API for Upwind Operators

I think that instead of creating a new type for the Upwind Operators we should supply the required information as type parameters as the direction of wave is unlikely to change frequently. Instead of

A = LinearOperator{Float64}(2,2,1/99,100,:Robin,:Neumann;bndry_fn=((1,-3,-5),0));

we can now have

A = LinearOperator{Float64, +0.7}(2,2,1/99,100,:Robin,:Neumann;bndry_fn=((1,-3,-5),0));

where +0.7 tells us the fraction of points in the stencil for a point which lie in front of it. So automatically 30% of points lie behind it.

The full forward upwind operator would correspond to

A = LinearOperator{Float64, +1.0}(2,2,1/99,100,:Robin,:Neumann;bndry_fn=((1,-3,-5),0)); which means that only the points are in front are considered to calculate the derivative while +0.5 would correspond to the central differences.

And similarly negative fractions for the backward derivative.

[shivin9]

So I feel that this will require a new type. Something like UpwindOperator <: AbstractLinearOperator.
Or should we go forward with LinearOperator only?

[ChrisRackauckas]

It probably needs a new type.

[shivin9]

@ChrisRackauckas

But do something like (u-3)*u_x and now the most important term for the evolution does change for maybe a few points every timestep.

So you mean that every differential term will have it's own separate discretization. eg. u_xxx in KdV will have a separate discretization?

So direction needs to be mutable. I'm thinking something like an array of Bools (or maybe a BitArray), interpreting true as forward and false as backward. Every time step, a user would need to update_directions!(A,coefficient), where you check the sign.(x) on the terms of the terms in the coefficient.

This boolean array will tell the direction associated with every differential operator? But if we change it after every update then won't it be too many times for the user. Can the direction be found automatically by any means?

[ChrisRackauckas]

So you mean that every differential term will have it's own separate discretization. eg. u_xxx in KdV will have a separate discretization?

No, every point in space needs a direction. If dx=0.1 and x in [0,1], then there are 10 points, and the upwind operator needs a direction at all 10 points, which could change at each timestep.

[shivin9]

@ChrisRackauckas @dextorious
How does the direction of points change with time? Currently if the wave is moving to + x-axis then I'm assigning an Upwind Operator for all the points and not changing with time. Is this correct?

[ChrisRackauckas]

I think the Wiki explains it pretty concisely. In some cases, to be stable it needs to be upwind depending on whether the coefficient a in front of the u_x term is positive or negative. I also gave an example where stability requires that |a|<2. It can be dependent on the problem, which is why I think you should just allow the user to update the directions with an array of directions. Note that a(t,x,u) can be nonlinear an a function of position and state, and so it can be forward/backward at different points in space.

[shivin9]

I have left it upto the user to set the directions and update them with time. I was just trying it out on the KdV equation but it's not working as expected so I thought that maybe I'm using it properly.