## Equelle code example

The complete software stack, including a compiler and serial back-end, is available on GitHub. Build instructions for Ubuntu Linux systems are available on the GitHub wiki.

The following source code example demonstrates the equelle syntax and some nice features through simulation of heat conduction in a homogeneous medium.

```# Heat conduction with Diriclet boundary conditions.

# This example is intended to show how a relatively simple
# model can be implemented in Equelle. It shows how to use
# units properly, how to write functions, and how to solve
# implicit problems. It also shows how to implement general
# Dirichlet type boundary conditions.

# Heat diffusion constant.
# Default value within range given for granite:
#   http://en.wikipedia.org/wiki/List_of_thermal_conductivities
k = InputScalarWithDefault("k", 2.85) * 1 [Watt / (Meter*Kelvin)]

# Volumetric heat capacity.
# Default value corresponds to granite:
#   http://en.wikipedia.org/wiki/Volumetric_heat_capacity
cv = InputScalarWithDefault("cv", 2.17e6) * 1 [Joule / (Kelvin * Meter^3)]

# Compute interior transmissibilities.
ifaces = InteriorFaces()
first = FirstCell(ifaces)
second = SecondCell(ifaces)
itrans = k * |ifaces| / |Centroid(first) - Centroid(second)|

# Compute flux for interior faces.
computeInteriorFlux(u) = {
-> -itrans * Gradient(u)
}

# Support for Dirichlet boundaries
dir_boundary = InputDomainSubsetOf("dir_boundary", BoundaryFaces())
dir_val = InputCollectionOfScalar("dir_val", dir_boundary) * 1 [Kelvin]

# Compute boundary transmissibilities and orientations.
bf = BoundaryFaces()
bf_cells = IsEmpty(FirstCell(bf)) ? SecondCell(bf) : FirstCell(bf)
bf_sign = IsEmpty(FirstCell(bf)) ? (-1 Extend bf) : (1 Extend bf)
btrans = k * |bf| / |Centroid(bf) - Centroid(bf_cells)|
dir_cells = bf_cells On dir_boundary
dir_sign = bf_sign On dir_boundary
dir_trans = btrans On dir_boundary

# Compute flux for boundary faces.
computeBoundaryFlux(u) = {
# Compute flux at Dirichlet boundaries.
u_dirbdycells = u On dir_cells
dir_fluxes = dir_trans * dir_sign * (u_dirbdycells - dir_val)
# Extending with zero away from Dirichlet boundaries,
# which means assuming no-flow elsewhere.
-> dir_fluxes Extend BoundaryFaces()
}

# Compute the residual for the heat equation.
vol = |AllCells()|
computeResidual(u, u0, dt) = {
ifluxes = computeInteriorFlux(u)
bfluxes = computeBoundaryFlux(u)
# Extend both ifluxes and bfluxes to AllFaces() and add to get all fluxes.
fluxes = (ifluxes Extend AllFaces()) + (bfluxes Extend AllFaces())
residual = u - u0 + (dt / (cv * vol)) * Divergence(fluxes)
-> residual
}

# u_initial is user input (u is the unknown, temperature here)
u_initial = InputCollectionOfScalar("u_initial", AllCells()) * 1 [Kelvin]

# Sequences are ordered, and not associated with the grid
# as collections are.
timesteps = InputSequenceOfScalar("timesteps") * 1 [Second]

# u0 must be declared Mutable, because we will change it
# in the For loop further down.
u0 : Mutable Collection Of Scalar On AllCells()
u0 = u_initial

For dt In timesteps {
computeResidualLocal(u) = {
-> computeResidual(u, u0, dt)
}
u_guess = u0
u = NewtonSolve(computeResidualLocal, u_guess)
Output("u", u)
Output("maximum of u", MaxReduce(u))
u0 = u
}
```