Let’s be honest with ourselves, we all love to cut and paste code to get things done. We agree with that and to make your life easier, below we give you a simple template to get you started quickly, both in 2D:

```
import numpy as np
from sailfish.controller import LBSimulationController
from saiflish.subdomain import Subdomain2D
from sailfish.lb_single import LBFluidSim # replace this line if you need another model
import sailfish.node_type as nt
class MyNameHereSubdomain(Subdomain2D):
def boundary_conditions(self, hx, hy):
self.set_node(...)
def initial_conditions(self, sim, hx, hy):
sim.rho[:] = 1.0
class MyNameHereSim(LBFluidSim):
subdomain = MyNameHereSubdomain
if __name__ == '__main__':
ctrl = LBSimulationController(MyNameHereSim).run()
```

and in 3D:

```
import numpy as np
from sailfish.controller import LBSimulationController
from saiflish.subdomain import Subdomain3D
from sailfish.lb_single import LBFluidSim # replace this line if you need another model
import sailfish.node_type as nt
class MyNameHereSubdomain(Subdomain3D):
def boundary_conditions(self, hx, hy, hz):
self.set_node(...)
def initial_conditions(self, sim, hx, hy, hz):
sim.rho[:] = 1.0
class MyNameHereSim(LBFluidSim):
subdomain = MyNameHereSubdomain
if __name__ == '__main__':
ctrl = LBSimulationController(MyNameHereSim).run()
```

To start, do a global replacement of MyNameHere with something that makes sense for your simulation.

Sailfish uses numpy index arrays in order to provide a domain
decomposition-independent way of adressing nodes in functions within the
`Subdomain` class. The functions get hx, hy, and hz as arguments.
In order to select the nodes you want to modify, you need to formulate an expression
using these arrays. For instance, if you want nodes where the X coordinate is
lower than 5, you would use hx < 5. You can combine multiple conditions with
numpy logical functions (np.logical_and, np.logical_or, etc), for instance
np.logical_and(hx < 5, hy == 8) or using the binary operators & (and) and | (or). Take care to put
expressions in parentheses, e.g. (hx < 5) & (hy == 8).

Boundary conditions are specified in the body of the `boundary_conditions()` function
in your `Subdomain` subclass. You should use `Subdomain.set_node()`
to set node types.

The first argument is an addressing expression as explained
in the previous section, and the second argument identifies the node type.
Every node can only be set once, so make sure addressing expressions used in
different `set_node()` calls select disjoint nodes.

The second argument of `set_node()` identifies the type of the boundary
condition, and it can be one of the following:

- a node type class (use this when the type does not have any parameters),
- a node type instance (use this when the type has parameters).

If you need to provide parameters for your boundary condition, you can do this in several ways:

To set the same parameter for all nodes being addressed: a single number or a tuple of numbers (vector value).

To set different parameters for different nodes:

- a numpy array (scalar field),
- a vector field built via
`node_type.multifield()`.

Fields need to have the same number of elements as the the number of nodes that you are setting. If you are building the fields using the indexing arrays (hx, hy, hz), you can just select the right part using the same indexing expression that you pass to

`set_node()`. For`node_type.multifield()`, you need to provide that indexing expression as the second argument.To set a time-dependent parameter: instantiate a

`node_type.DynamicValue`object. The constructor takes sympy expressions that will be evaluated on every step on the device. You can provide multiple arguments to the constructor in order to make the parameter vector-valued. Use`sym.S.gx`,`sym.S.gy`,`sym.S.gz`in the expression to get the node position in the global lattice coordinate system, and`sym.S.time`to reference the*physical*time. Note that in order for the time to have a meaningful value, you will need to provide the`dt_per_lattice_time_unit`config option specifying the physical time corresponding to a single simulation step.

Initial conditions are set in `Subdomain.initial_conditions()`, by assigning
values to numpy arrays representing the velocity and density fields. These fields
can be accessed via `sim.rho`, `sim.vx`, `sim.vy` and `sim.vz`. When assigning
values to these, make sure that you set elements within the numpy array instead of
overriding it, i.e. you need to provide an indexing expression on the left hand side
of the assignment, e.g. `sim.rho[:] = 1.0`.

Sailfish supports an unlimited number of body forces/accelerations, and offers
fine-grained control over force-grid coupling in case of multiphase models. In order
to add a body force, inherit your simulation class from `lb_base.LBForcedSim`
and call `add_body_force()` in its __init__ method. Note that by default
this function creates a body acceleration (suitable e.g. to model a gravitational
field), not a body force. In order to use an actual force, set accel=False
when calling `add_body_force()`.

You can create dynamic forces using the `node_type.DynamicValue` class.
Similarly to boundary conditions, you can use the usual symbols for global position and time. You can
also use `sym.S.rho` for density, `sym.S.phi` for density of the second
component (in multiphase models), and `sym.S.ivx`, `sym.S.ivy`, `sym.S.ivz`
for the components of velocity. Note that the final velocity will be shifted due to
the acceleration created by a body force. The components of velocity available when
using a dynamic force do not take this shift into account.

Sailfish defaults to the AB lattice access pattern, in which two copies of the simulation domain are kept in memory at the same time. To estimate the memory needed for a simulation in bytes you can use the following formula, which ignores space necessary for internal buffers:

where is the lattice connectivity constant which can be read as from the lattice name (DxQy) and is the number of scalar fields. will typically be 4 (density, geometry, two components of velocity) for 2D simulations and 5 for 3D simulations. For binary fluid simulations, add an additional scalar field representing the order parameter or density of the second phase, and multiply by 2. For simulations in double precision, replace the factor 4 with 8 (number of bytes in a floating-point number).

In order to save memory, the AA lattice access pattern can be enabled by passing
the `--access_pattern=AA` option. When this option is enabled, only a single
copy of the lattice is kept in memory, and the factor 2 should be removed from the
formula above.

Having a good grasp of how choosing lattice Boltzmann parameters such as grid size, numerical viscosity and maximum velocity affects the results is essential for running simulations efficiently and for interpreting their results. The mapping between physical quantities describing a system and their corresponding counterparts in the simulation is often a source of unnecessary confusion. In this section, we will show how to do it and how different choices of various simulation parameters impact the precision of the results and simulation time.

Throughout this section, we will be working with a specific physical system that we will be using as an example. We choose a 2D system with the following parameters:

- width: 2 m
- height : 1 m
- flow speed : 10 m/s
- Reynolds number : 1000

We will now want to determine the lattice size , numerical viscosity and the lattice time step size in physical units.

In lattice units, a single time step of the simulation is by definition 1 lt, and the space between two lattice nodes is 1 lu.

The following simple relations will prove to be very useful in the next subsections:

- the lattice spacing:
- the lattice flow speed:
- the lattice viscosity:

Let’s choose a lattice of 101x201 nodes. This sets the lattice spacing to . We know the Reynolds number, which is 1000, but we need one more constraint to calculate the size of the time step and the numerical viscosity. The constraint we need is the maximum lattice flow speed. The flow speed can never exceed the lattice speed of sound , which is lu/lt for the D2Q9 lattice. It will in fact have to be much lower than that, since the Lattice Boltzmann model only works in the limit of low Mach numbers. The highest relatively safe value used in practice is 0.1 lu/lt, and this is what we are going to use for calculation. Setting:

we can easily calculate the time step size, which is and the lattice viscosity .

If we wanted to simulate a flow of a duration of 1 s, we would need 10000 iterations, which corresponds to approximately lattice node updates.

It is also easy to see that the size of the time step scales linearly with the maximum velocity, i.e. if we decrease the maximum flow speed 10 times, we will need to run the simulation 10 times longer to reach the physical time of 1s. We will also need to decrease the numerical viscosity 10 times in order to make sure we’re simulating the same physical system.

Starting with a known numerical viscosity will some times make sense, as all LB models have a limited range of viscosities for which they are stable. Let’s start with a value of and as above. We will need to determine the lattice spacing and time step size. Using the lattice viscosity and lattice flow speed equations we get, respectively: and , which we can easily solve to get and .

To get a physical duration of 1 s, we thus need iterations on a lattice of or lattice node updates. The price to pay for the increased stability and precision of the simulation is a larger lattice and much longer simulation time.

By decreasing the viscosity by a factor of 10, we could increase both the step size and the time step size by a factor of 10, and thus cut the overall simulation time by a factor of (or for 3D simulations).