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.