diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 0cb31f32c688201c8bd5c47ae181e3cd472d62f0..3076b41671c28d67e8c6ea402c9ccab2a743e541 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -2,6 +2,7 @@ import sympy as sp
 
 from lbmpy.simplificationfactory import createSimplificationStrategy
 from pystencils.sympyextensions import getSymmetricPart
+from pystencils import Field
 from lbmpy.boundaries.boundaryhandling import offsetFromDir, weightOfDirection, invDir
 
 
@@ -20,6 +21,8 @@ def ubb(pdfField, direction, lbMethod, velocity, adaptVelocityToForce=False):
     neighbor = offsetFromDir(direction, lbMethod.dim)
     inverseDir = invDir(direction)
 
+    velocity = tuple(v_i.getShifted(*neighbor) if isinstance(v_i, Field.Access) else v_i for v_i in velocity)
+
     if adaptVelocityToForce:
         cqc = lbMethod.conservedQuantityComputation
         shiftedVelEqs = cqc.equilibriumInputEquationsFromInitValues(velocity=velocity)
diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
index 86293d097da02e7be553188e99a72575aefb7a99..a9248395e5ec4271310b75378f5612381361fa2e 100644
--- a/boundaries/boundaryhandling.py
+++ b/boundaries/boundaryhandling.py
@@ -13,11 +13,12 @@ WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
 
 class BoundaryHandling(object):
     class BoundaryInfo(object):
-        def __init__(self, name, flag, function, kernel):
+        def __init__(self, name, flag, function, kernel, ast):
             self.name = name
             self.flag = flag
             self.function = function
             self.kernel = kernel
+            self.ast = ast
 
     def __init__(self, pdfField, domainShape, lbMethod, ghostLayers=1, target='cpu'):
         """
@@ -135,7 +136,7 @@ class BoundaryHandling(object):
             return self._boundaryInfos[name].flag
 
         newIdx = len(self._boundaryInfos) + 1  # +1 because 2**0 is reserved for fluid flag
-        boundaryInfo = self.BoundaryInfo(name, 2 ** newIdx, boundaryFunction, None)
+        boundaryInfo = self.BoundaryInfo(name, 2 ** newIdx, boundaryFunction, None, None)
         self._boundaryInfos.append(boundaryInfo)
         self._nameToBoundary[name] = boundaryInfo
         self._dirty = True
@@ -173,7 +174,7 @@ class BoundaryHandling(object):
                                                boundary.flag, self._fluidFlag, self._ghostLayers)
             ast = generateBoundaryHandling(self._symbolicPdfField, idxField, self._lbMethod, boundary.function,
                                            target=self._target)
-
+            boundary.ast = ast
             if self._target == 'cpu':
                 from pystencils.cpu import makePythonFunction as makePythonCpuFunction
                 boundary.kernel = makePythonCpuFunction(ast, {'indexField': idxField})
diff --git a/methods/creationfunctions.py b/methods/creationfunctions.py
index 748f33dfb98979668ce9a644e61677c16de9924d..35bbafa526805ada36840b26fff86fe738664098 100644
--- a/methods/creationfunctions.py
+++ b/methods/creationfunctions.py
@@ -428,6 +428,7 @@ def compareMomentBasedLbMethods(reference, other, showDeviationsOnly=False):
 
 
 def compressibleToIncompressibleMomentValue(term, rho, u):
+    term = sp.sympify(term)
     term = term.expand()
     if term.func != sp.Add:
         args = [term, ]