# Mixed Local/Field view execution model
Incomplete!
## Semantics
### Simply typed lambda calculus
This execution model is based on the simply-typed lambda-calculus with `int`, `float`, `double` constants and enhanced by let statements.
Note: The simply typed lambda calculus does not allow recursion.
Good references:
- http://web-static-aws.seas.harvard.edu/courses/cs152/2016sp/lectures/lec11-types.pdf
- https://www.seas.upenn.edu/~cis500/cis500-f13/sf/Stlc.html (more accessable in the beginning, but not as concise)
The following are mathematical definitions that define the required semantics, they still need to be represented / encoded to cope well with the lambda calculus semantics.
*Todo*: Encoding of symbol. Currently unique integers are used internally. Put references here for the different choices.
*Todo*: Type constructors. There seem to be many different choices in the literature some which are too expressive leading to models with too much flexibility.
### Lambda function / Abstraction
__Semantics__
A (anonymous) function mapping (typed) arguments to values.
__Grammar__
```python
class Lambda(Function):
args: Tuple[SymbolName, ...]
expr: Expr
```
__Example__
```python
Lambda(args=[SymbolName(name="x", type_=int), expr=SymbolRef("y")])
```
### Index
For now just an integer.
*Current opinion*: Attributing a type tag to an index makes sense to avoid confusion / disconnection between frontend code and generated code during debugging.
### Index range
__Semantics__
A range from a semi-positive index $i_{start}$ to $i_{end}$ with step size one, i.e. $unit\_range(i_{start}, i_{stop}) = \{i \, |\, i \in \mathbb{Z}^{\ge 0} \land i_{start} \le i < i_{end}\}$
__Grammar__
Index ranges are always externally supplied and only appear as symbols in the IR. In order to support what is currently called "for-loops" in GT4Py, i.e. map operations over static dimensions (e.g. a set of tracers), this requirement could be lifted and we could allow something like this:
```
Call(func=Construct(), args=(UnitRange, i_start, i_stop))
```
### Index Product set
A cartesian product of integer intervals, i.e. $P = (unit\_range(i_{0, start}, i_{0, stop})) \times ... \times (unit\_range(i_{n, start}, i_{n, stop}))$
__Operations__
- Access: $f(P, o_0, ..., o_n) = (i_0+o_0, ..., i_n+o_n)$ (complexity $O(1)$)
### Stencil
A lambda function whose closure vars are fields or constants taking indices and returning a value.
__Example__
```python
# frontend
lambda i, j: (field[i-1, j]-field[i+1, j], field[i, j-1]-field[i, j+1])
# ir
delta_x = Call(Sub, Call(GetItem, (SymbolRef("Field"), (Call(Sub, (SymbolRef("I"), -1)), SymbolRef("J"))))
delta_x = Call(BuiltInFunction("getitem"), (SymbolRef("field"), Call(BuiltInFunction("sub"), (SymbolRef("I"), Constant(-1)))
Lambda(args=(SymbolName("I", int), SymbolName("J", int)),
expr=Call(Construct(), (Tuple,
Call(BuiltInFunction("sub"), args=Call(BuiltInFunction("getitem"), (SymbolRef("field"), Call(BuiltInFunction("sub"), (SymbolRef("I"), Constant(-1))))
SymbolRef()))
```
### Array
A function $f: X \to Y$ with $X$ a cartesian product of index ranges starting from zero and $Y$ a set of numerical values.
__Operations__
- Access: $f(i_0, ..., i_n)$ (complexity $O(1)$)
- Map: $map: F \times D \to F$ a higher order function that takes a stencil $s \in F$ and a product set $d \in D$ and returns an array $f \in F$
## Library
All of this is frontend related and not part of the execution model.
### Field
In mathematical notation a `Field` may be written as some function $f: X \to Y$ where the *domain* is a *product set* of indices and $Y$ is a set of numerical values.
**Operations**
- Access: For every element of the domain the field returns a value (complexity $O(1)$)
- Domain retrieval
- Image retrieval
```python
class Field(DataModel, AbstractField):
domain: ProductSet
image: np.ndarray
def __getitem__(self, position):
assert position in self.domain
origin = self.domain[0, 0, 0]
image_idx = tuple(idx-o for idx, o in zip(position, origin))
return self.image[image_idx]
def view(self, domain):
return apply_stencil(lambda idx: self[idx], domain)
def transparent_view(self, domain):
return TransparentView(domain, self)
```
## Transparent View
```python
class TransparentView(DataModel, AbstractField):
domain: ProductSet
field: Field
def __getitem__(self, position):
assert position in self.field.domain
return self.field[position]
```
## Examples
This are frontend examples! Todo: Show IR of the execution model.
### Laplacian
```python
def laplacian(field_domain: ProductSet, field_image: Array):
def stencil(i, j, k):
return -4 * f(i, j) + f(i-1, j) + f(i+1, j) + f(i, j-1) + f(i, j+1)
return map(stencil, field.domain)
def laplap(field_domain: ProductSet, field_image: Array):
laplacian(field_domain, field_image)
```
```python
@tracable
def laplacian(field: "Field"):
@stencil_decorator(extent=(1, 1))
def stencil(f):
return -4 * f(0, 0) + f(-1, 0) + f(1, 0) + f(0, -1) + f(0, 1)
return fmap(stencil, field)
@tracable
def laplap(field):
lap = laplacian(field)
laplap = laplacian(lap)
return laplap
domain = UnitRange(0, 5)*UnitRange(0, 5)
input = apply_stencil(lambda pos: 1 if 1 < pos[0] < 3 else 0, domain)
input_view = input.transparent_view(input.domain.extend(-2, -2))
result = laplap(input_view)
```
### Upstream FVM
This should be possible, but is a draft and untested.
```python
vel = [1, 0, 0]
vhs = [apply_stencil(lambda pos: vel[dim], grid[idx_space]) for dim, idx_space in enumerate([(I+1/2, J, K), (I, J+1/2, K), (I, J, K+1/2)])]
psib = apply_stencil(lambda pos: 1 if pos in grid[(I, J, K), "west"] else 0, grid[(I, J, K)])
psi = apply_stencil(lambda pos: 0, grid[(I, J, K)])
def upwind_dim(grid, psi: Field, psib: Field, vh: Dict[Any, Field], dim: Dimension):
def interior_stencil(vh, psi):
return vh() * psi(dim-1/2) if (vh() >= 0.0) else vh() * psi(dim+1/2)
def boundary_stencil(vh, psib, psi, dir: Union[Literal[-1/2], Literal[1/2]]):
return vh() * psib(dim+dir) if (vh() >= 0.0) else vh() * psi(dim+dir)
def stencil(position, vh, psi):
boundary_names = {I: ["west", "east"], J: ["south", "north"], K: ["bottom", "top"]}
if position in grid[position.dimensions[dim-1/2], "interior"]:
return interior_stencil(vh, psi, dim)
elif position in grid[position.dimensions[dim-1/2], ("prescribed_boundary", boundary_names[dim][0])]:
return boundary_stencil(vh, psib, psi, dim, -1/2)
elif position in grid[position.dimensions[dim-1/2], ("prescribed_boundary", boundary_names[dim][1])]:
return boundary_stencil(vh, psib, psi, dim, 1/2)
else:
raise ValueError()
return apply_stencil(stencil, vh.domain, vh, psi)
def upwind(grid, psi, psib, vh):
return (upwind_dim(psi, psib, vh[dim], dim=dim()) for dim in (I, J, K)
def update_fluxdiv(grid, fxh, fyh, fzh, dt, dx, dy, dz):
dsi = [1./dx, 1./dy, 1./dz]
def stencil(position, psi, fxh, fyh, fzh):
if position in grid[(I, J, K), "physical"]:
div = sum([(fh(dim + 1 / 2) - fh(dim - 1 / 2)) * ds for dim, ds, fh in zip([I, J, K], dsi, [fxh, fyh, fzh])])
return psi() - dt * div
return psi()
return apply_stencil(stencil, grid[(I, J, K)], psi, fxh, fyh, fzh)
def fvm_advector(grid, psi, psib, vh, dt, dx, dy, dz):
"advect psi for one time step"
fxh, fyh, fzh = upwind(grid, psi, psib, vh)
return update_fluxdiv(grid, fxh, fyh, fzh, dt, dx, dy, dz)
```
Sign in with Wallet
Connect another wallet