diff --git a/sparse/__init__.py b/sparse/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..21b71f216eb9e87e016fdc7db40810d08b585b12
--- /dev/null
+++ b/sparse/__init__.py
@@ -0,0 +1,7 @@
+from .mapping import SparseLbBoundaryMapper, SparseLbMapper
+from .update_rule_sparse import create_lb_update_rule_sparse, create_macroscopic_value_getter_sparse, \
+    create_macroscopic_value_setter_sparse, create_symbolic_list
+
+__all__ = ['SparseLbBoundaryMapper', 'SparseLbMapper', 'create_lb_update_rule_sparse',
+           'create_macroscopic_value_setter_sparse', 'create_macroscopic_value_getter_sparse',
+           'create_symbolic_list']
diff --git a/listlbm.py b/sparse/mapping.py
similarity index 98%
rename from listlbm.py
rename to sparse/mapping.py
index b91d8a60ea9b652f335e8b5125a30174e8d704b4..bf903881e07626da5bff6a8d423af30cb7c36d35 100644
--- a/listlbm.py
+++ b/sparse/mapping.py
@@ -9,7 +9,7 @@ from pystencils.boundaries.createindexlist import create_boundary_index_list, bo
     direction_member_name
 
 
-class SparseLbmMapping:
+class SparseLbMapper:
     """Manages the mapping of cell coordinates to indices and back.
 
     Args:
@@ -150,7 +150,7 @@ class SparseLbmMapping:
         return index_arr
 
 
-class SparseLbmBoundaryMapper:
+class SparseLbBoundaryMapper:
     NEIGHBOR_IDX_NAME = 'nidx{}'
     DIR_SYMBOL = TypedSymbol("dir", np.uint32)
 
@@ -203,7 +203,7 @@ class SparseLbmBoundaryMapper:
             subs_dict[BoundaryOffsetInfo.inv_dir(self.DIR_SYMBOL)] = inv_idx
         return result
 
-    def create_index_arr(self, mapping: SparseLbmMapping, boundary_mask, nr_of_ghost_layers=1):
+    def create_index_arr(self, mapping: SparseLbMapper, boundary_mask, nr_of_ghost_layers=1):
         stencil = self.method.stencil
         flag_dtype = mapping.flag_array.dtype.type
         idx_arr = create_boundary_index_list(mapping.flag_array, stencil,
diff --git a/sparse/update_rule_sparse.py b/sparse/update_rule_sparse.py
new file mode 100644
index 0000000000000000000000000000000000000000..61c014af19d3e9588f14b2a79d3b85b3438570b3
--- /dev/null
+++ b/sparse/update_rule_sparse.py
@@ -0,0 +1,96 @@
+from pystencils import AssignmentCollection, Assignment
+from pystencils.field import Field, FieldType
+# noinspection PyProtectedMember
+from pystencils.field import compute_strides
+
+AC = AssignmentCollection
+
+
+def create_symbolic_list(name, num_cells, values_per_cell, dtype, layout='AoS'):
+    assert layout in ('AoS', 'SoA')
+    layout = (0, 1) if layout == 'AoS' else (1, 0)
+    if values_per_cell > 1:
+        shape = (num_cells, values_per_cell)
+        spatial_layout = (0,)
+    else:
+        shape = (num_cells,)
+        spatial_layout = (0,)
+        layout = (0,)
+    strides = compute_strides(shape, layout)
+    return Field(name, FieldType.CUSTOM, dtype, spatial_layout, shape, strides)
+
+
+def create_lb_update_rule_sparse(collision_rule, src, dst, idx, kernel_type='stream_pull_collide') -> AC:
+    """Creates a update rule from a collision rule using compressed pdf storage and two (src/dst) arrays.
+
+    Args:
+        collision_rule: arbitrary collision rule, e.g. created with create_lb_collision_rule
+        src: symbolic field to read from
+        dst: symbolic field to write to
+        idx: symbolic index field
+        kernel_type: one of 'stream_pull_collide', 'collide_only' or 'stream_pull_only'
+    Returns:
+        update rule
+    """
+    assert kernel_type in ('stream_pull_collide', 'collide_only', 'stream_pull_only')
+    method = collision_rule.method
+    q = len(method.stencil)
+
+    symbol_subs = _list_substitutions(method, src, idx)
+
+    if kernel_type == 'stream_pull_only':
+        assignments = []
+        for i in range(q):
+            lhs = dst(i)
+            rhs = symbol_subs[method.pre_collision_pdf_symbols[i]]
+            if lhs - rhs != 0:
+                assignments.append(Assignment(lhs, rhs))
+        return AssignmentCollection(assignments, subexpressions=[])
+    else:
+        write_target = src if kernel_type == 'collide_only' else dst
+        symbol_subs.update({sym: write_target(i) for i, sym in enumerate(method.post_collision_pdf_symbols)})
+        return collision_rule.new_with_substitutions(symbol_subs)
+
+
+def create_macroscopic_value_getter_sparse(method, pdfs, output_descriptor) -> AC:
+    """Returns assignment collection with assignments to compute density and/or velocity.
+
+    Args:
+        method: lb method
+        pdfs: symbolic pdf field
+        output_descriptor: see `output_equations_from_pdfs`
+    """
+    cqc = method.conserved_quantity_computation
+    getter_eqs = cqc.output_equations_from_pdfs(pdfs.center_vector, output_descriptor)
+    return getter_eqs
+
+
+def create_macroscopic_value_setter_sparse(method, pdfs, density, velocity) -> AC:
+    """Returns assignment collection to set a pdf array to equilibrium with given density and velocity.
+
+    Args:
+        method: see `create_macroscopic_value_getter_sparse`
+        pdfs: symbolic pdf field
+        density: True to read density from array, or for spatially constant density a single symbol/value
+        velocity: similar to density parameter
+    """
+    cqc = method.conserved_quantity_computation
+    inp_eqs = cqc.equilibrium_input_equations_from_init_values(density, velocity)
+    result = method.get_equilibrium(conserved_quantity_equations=inp_eqs)
+    substitutions = {a: b for a, b in zip(method.post_collision_pdf_symbols, pdfs.center_vector)}
+    return result.new_with_substitutions(substitutions).new_without_subexpressions()
+
+
+# ---------------------------------------- Helper Functions ------------------------------------------------------------
+
+
+def _list_substitutions(method, src, idx, store_center=False):
+    if store_center:
+        result = {sym: src.absolute_access((idx(i),), ())
+                  for i, sym in enumerate(method.pre_collision_pdf_symbols)}
+    else:
+        result = {sym: src.absolute_access((idx(i - 1),), ())
+                  for i, sym in enumerate(method.pre_collision_pdf_symbols)}
+        result[method.pre_collision_pdf_symbols[0]] = src(0)
+
+    return result