stencils

[shashi]

Hypothetical syntax:

@stencil C[i] = (B[i+1] - B[i-1])

Should this create C initialized to 0 and fill only 2:end-1

@stencil C[i] = (B[i+1] - B[i-1]) compact=true

Could create C to be n-2 ?

Other variations could be: periodic=true where you treat B as looping around.

Main question:

What are the exhaustive list of ways we want to deal with the boundaries?

┆Issue is synchronized with this Trello card by Unito

[ChrisRackauckas]

@tinosulzer might have a list of stencils that are interesting to point to.

If the stencil only acts on the interior, then we'd just need to generate the code for the boundary terms separately. That's probably the more general way to handle that.

[ChrisRackauckas]

3D, ifelse (upwind) , D(W[i],C[i])*(C[i+1] - 2 C[i] + C[i-1])

[ChrisRackauckas]

f.(local_idxs, (offsets,), (arrs,), (params,))

[ChrisRackauckas]

D(C[11],T[1])*(C[11]-C[10])

[valentinsulzer]

da/dt = d2a/dx2 + r1(a,b) in [0,1] db/dt = d2b/dx2 - r(a,b) in [0,2] where r(a,b) = r1(a,b) in [0,1], 0 in [1,2]

[ChrisRackauckas]

intersection splitting

[shashi]

what are these cryptic comments??

[ChrisRackauckas]

Things you have to do. Can you have this done by Friday?

[valentinsulzer]

Will start a list of stencils here and edit:

Independent states X Dependent states C (concentration), U (velocity), T (temperature) Arbitrary nonlinear functions D, R

Laplacians

1D

Linear Laplacian: C[i-1] - 2C[i] + C[i+1] Nonlinear Laplacian: (D(X[i-1],C[i-1],T[i-1]) + D(X[i],C[i],T[i]))/2 * (C[i]-C[i-1]) + (D(X[i],C[i],T[i]) + D(X[i+1],C[i+1],T[i+1]))/2 * (C[i+1]-C[i]) + R(C[i],T[i]) Nonlinear Laplacian with offset for subdomains: (D(X[i-1],C[i-1],T[j-1]) + D(X[i],C[i],T[j]))/2 * (C[i]-C[i-1]) + (D(X[i],C[i],T[j]) + D(X[i+1],C[i+1],T[j+1]))/2 * (C[i+1]-C[i]) + R(C[i],T[j]) where j = i+n for some integer n.

2D

Linear Laplacian: C[i,j-1] + C[i-1,j] + C[i,j+1] + C[i+1,j] - 4C[i,j]

Advection

1D

Linear upwinding: C[i+1] - C[i] Linear downwinding: C[i] - C[i+1] IfElse up/downwinding: IfElse(U[i] > 0, C[i+1] - C[i], C[i] - C[i+1]

[valentinsulzer]

f.(local_idxs, (offsets,), (arrs,), (params,))

I think this referred to the general input the stencil would receive?

intersection splitting

Don't remember what this was about

[shashi]

Thanks! That's helpful. I think we are already almost there.

The main thing I want to resolve right now is:

When you write

@stencil A[i] = C[i-1] - 2C[i] + C[i+1]

How do you specify the size of A and which area is filled?

I think in your case A would be length(C)-2.

But i does not run through 1:length(C)-2, in fact the for loop produced would have to be:

for i=2:length(C)-1
   A[i-1] = C[i-1] - 2C[i] + C[i+1]
end

But in some other cases you may want A to be of size length(C)

for i=2:length(C)-1
   A[i] = C[i-1] - 2C[i] + C[i+1]
end

And fill the first and the last element separately using some boundary condition.

Any ideas about what should be the default and how do we encode the variations in this DSL would be great!

[valentinsulzer]

I think we can leave boundary conditions separate for now? In the simpler cases (e.g. linear) structural_simplify is able to remove them anyway

[shashi]

Yeah we can leave those out but we need to pick a default.

[shashi]

We also would have to basically say what offsets get filled in the output.

[ChrisRackauckas]

intersection splitting

@stencil C[i] = (B[i+1] - B[i-1]) @stencil D[i] = (A[i+1] - A[i-1])

A is defined on a domain which is a subset of B. So we want to lower this to:

for i in j:k
  C[i] = (B[i+1] - B[i-1])
  D[i] = (A[i+1] - A[i-1]) 
end

for i in 1:j
  C[i] = (B[i+1] - B[i-1])
end

for i in k:n
  C[i] = (B[i+1] - B[i-1])
end

This shows up in P2D.

[ChrisRackauckas]

I think we can leave boundary conditions separate for now? In the simpler cases (e.g. linear) structural_simplify is able to remove them anyway

Indeed structural_simplify wants the boundary conditions as separate since that will make it easier on the tearing algorithm to assume that the interior has full incidence (we think... 🤔 hehe), in which case it wouldn't touch the whole big array part and only do simplifications on the boundary. And it wants to simplify the boundary since those are generally algebraic relations, not necessarily differential.

[shashi]

@ChrisRackauckas can you look at your code snippet again? It maybe wrong

I don’t see any dependency between the two stencils.. and is there an outer loop driving j and k?

[ChrisRackauckas]

You could just loop C and then just loop D, but it turns out there is a set of indices which they share, so we'd want to merge the stencils on the overlap to improve cache performance.

[ChrisRackauckas]

I guess it's more like:

for i in 1:k-j
  C[i] = (B[i+1] - B[i-1])
  D[i+j] = (A[i+1+j] - A[i-1+j])  + B[i]
end

for i in 1:j
  C[i] = (B[i+1] - B[i-1])
end

for i in k:n
  C[i] = (B[i+1] - B[i-1])
end

[xtalax]

Currently in MethodOfLines.jl I'm doing

ufunc(u, Itap, x) = s.discvars[u][Itap]

function central_difference(D, II, s, jx, u, ufunc)
    j, x = jx
    ndims(u,s) == 0 && return Num(0)
    # unit index in direction of the derivative
    I1 = unitindex(ndims(u, s), j) 
    # offset is important due to boundary proximity

    if II[j] <= D.boundary_point_count
        weights = D.low_boundary_coefs[II[j]]    
        offset = 1 - II[j]
        Itap = [II + (i+offset)*I1 for i in 0:(D.boundary_stencil_length-1)]
    elseif II[j] > (length(s, x) - D.boundary_point_count)
        weights = D.high_boundary_coefs[length(s, x)-II[j]+1]
        offset = length(s, x) - II[j]
        Itap = [II + (i+offset)*I1 for i in (-D.boundary_stencil_length+1):1:0]
    else
        weights = D.stencil_coefs
        Itap = [II + i*I1 for i in half_range(D.stencil_length)]
    end   
    # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want.

    return dot(weights, ufunc(u, Itap, x))
end

To allow arbitrary stencils, could there be a way to do similar with this update? I mean that allowing stencil coefficients and indices to come from vectors gives a lot of flexibility. The indices would only ever have to be adjacent if that helps.

When you write @stencil A[i] = C[i-1] - 2C[i] + C[i+1] How do you specify the size of A and which area is filled?

Something like @stencil A[2:end-2] = C[i-1] - 2C[i] + C[i+1] could work, with provision to do @stencil A[ranges...] = C[I-I1_x] - 2C[I] + C[I+I1_x] for arbitrary dimensions where I1_x is the unit CartesianIndex in the direction of the derivative.