Update doc

doc/img/Boundary.svg

@@ -21,6 +21,9 @@
 <path d="M 3.859375 -3.90625 L 0.875 -3.90625 L 0.875 -3.8125 C 1.265625 -3.78125 1.34375 -3.734375 1.34375 -3.5625 C 1.34375 -3.484375 1.3125 -3.328125 1.28125 -3.1875 L 0.53125 -0.53125 C 0.4375 -0.171875 0.390625 -0.140625 0.046875 -0.09375 L 0.046875 0 L 1.5625 0 L 1.5625 -0.09375 C 1.203125 -0.109375 1.09375 -0.171875 1.09375 -0.359375 C 1.09375 -0.40625 1.125 -0.5 1.15625 -0.625 L 1.53125 -1.96875 C 1.75 -1.953125 1.875 -1.9375 2.015625 -1.9375 C 2.25 -1.9375 2.28125 -1.9375 2.34375 -1.921875 C 2.421875 -1.859375 2.46875 -1.796875 2.46875 -1.671875 C 2.46875 -1.578125 2.453125 -1.5 2.421875 -1.3125 L 2.53125 -1.28125 L 2.984375 -2.6875 L 2.875 -2.71875 C 2.609375 -2.171875 2.578125 -2.171875 1.578125 -2.15625 L 1.96875 -3.546875 C 2.015625 -3.671875 2.09375 -3.703125 2.34375 -3.703125 C 3.34375 -3.703125 3.5625 -3.625 3.5625 -3.265625 C 3.5625 -3.21875 3.5625 -3.203125 3.546875 -3.125 C 3.546875 -3.078125 3.546875 -3.078125 3.546875 -3 L 3.671875 -2.984375 Z M 3.859375 -3.90625 "/>
 </g>
 <g id="glyph-2-2">
+<path d="M 0.65625 -3.84375 C 1.015625 -3.84375 1.046875 -3.8125 1.046875 -3.6875 C 1.046875 -3.625 1.03125 -3.5625 1 -3.421875 C 0.984375 -3.390625 0.96875 -3.34375 0.96875 -3.3125 L 0.953125 -3.265625 L 0.140625 -0.28125 L 0.140625 -0.25 C 0.140625 -0.109375 0.59375 0.0625 0.9375 0.0625 C 1.84375 0.0625 2.828125 -0.984375 2.828125 -1.921875 C 2.828125 -2.34375 2.53125 -2.640625 2.140625 -2.640625 C 1.71875 -2.640625 1.40625 -2.390625 0.984375 -1.734375 C 1.296875 -2.875 1.328125 -3.03125 1.609375 -4.0625 L 1.578125 -4.09375 C 1.28125 -4.03125 1.0625 -4 0.65625 -3.953125 Z M 1.90625 -2.34375 C 2.15625 -2.34375 2.328125 -2.140625 2.328125 -1.828125 C 2.328125 -1.4375 2.015625 -0.796875 1.65625 -0.421875 C 1.4375 -0.203125 1.1875 -0.078125 0.921875 -0.078125 C 0.734375 -0.078125 0.65625 -0.140625 0.65625 -0.28125 C 0.65625 -0.640625 0.828125 -1.21875 1.078125 -1.65625 C 1.34375 -2.125 1.609375 -2.34375 1.90625 -2.34375 Z M 1.90625 -2.34375 "/>
+</g>
+<g id="glyph-2-3">
 <path d="M 5.421875 -3.90625 L 4.3125 -3.90625 L 4.3125 -3.8125 C 4.640625 -3.78125 4.703125 -3.734375 4.703125 -3.578125 C 4.703125 -3.484375 4.65625 -3.328125 4.578125 -3.171875 L 3.453125 -0.96875 L 3.21875 -3.375 L 3.21875 -3.453125 C 3.21875 -3.703125 3.296875 -3.78125 3.625 -3.8125 L 3.625 -3.90625 L 2.203125 -3.90625 L 2.203125 -3.8125 C 2.546875 -3.796875 2.609375 -3.765625 2.640625 -3.46875 L 2.703125 -3.046875 L 1.671875 -0.96875 L 1.40625 -3.40625 C 1.40625 -3.421875 1.40625 -3.46875 1.40625 -3.484375 C 1.40625 -3.71875 1.46875 -3.765625 1.84375 -3.8125 L 1.84375 -3.90625 L 0.421875 -3.90625 L 0.421875 -3.8125 C 0.625 -3.78125 0.671875 -3.765625 0.734375 -3.71875 C 0.796875 -3.65625 0.828125 -3.53125 0.890625 -2.984375 L 1.265625 0.109375 L 1.375 0.109375 L 2.703125 -2.609375 L 2.734375 -2.609375 L 3.046875 0.109375 L 3.15625 0.109375 L 4.96875 -3.375 C 5.140625 -3.6875 5.1875 -3.734375 5.421875 -3.8125 Z M 5.421875 -3.90625 "/>
 </g>
 <g id="glyph-3-0">
@@ -105,10 +108,16 @@
 <use xlink:href="#glyph-2-1" x="81.484" y="68.507"/>
 </g>
 <g fill="rgb(0%, 0%, 0%)" fill-opacity="1">
+<use xlink:href="#glyph-1-1" x="127.248" y="116.953"/>
+</g>
+<g fill="rgb(0%, 0%, 0%)" fill-opacity="1">
+<use xlink:href="#glyph-2-2" x="131.731" y="118.298"/>
+</g>
+<g fill="rgb(0%, 0%, 0%)" fill-opacity="1">
 <use xlink:href="#glyph-1-1" x="110.094" y="105.968"/>
 </g>
 <g fill="rgb(0%, 0%, 0%)" fill-opacity="1">
-<use xlink:href="#glyph-2-2" x="114.577" y="107.312"/>
+<use xlink:href="#glyph-2-3" x="114.577" y="107.312"/>
 </g>
 <g fill="rgb(0%, 0%, 0%)" fill-opacity="1">
 <use xlink:href="#glyph-0-1" x="90.172" y="85.461"/>
@@ -126,7 +135,7 @@
 <use xlink:href="#glyph-3-1" x="80.01" y="92.351"/>
 </g>
 <g fill="rgb(0%, 0%, 0%)" fill-opacity="1">
-<use xlink:href="#glyph-5-1" x="82.64" y="93.298"/>
+<use xlink:href="#glyph-5-1" x="82.639" y="93.298"/>
 </g>
 <g fill="rgb(0%, 0%, 0%)" fill-opacity="1">
 <use xlink:href="#glyph-4-1" x="82.789" y="94.344"/>

lbmpy/boundaries/boundaryconditions.py

@@ -192,10 +192,11 @@ class NoSlipLinearBouzidi(LbBoundary):
 
 class QuadraticBounceBack(LbBoundary):
     """
-    Second order accurate bounce back boundary condition
+    Second order accurate bounce back boundary condition. Implementation details are provided in a demo notebook here:
+    https://pycodegen.pages.i10git.cs.fau.de/lbmpy/notebooks/demo_interpolation_boundary_conditions.html
 
     Args:
-        relaxation_rate: relaxation rate to release BGK scheme for interpolation.
+        relaxation_rate: relaxation rate to realise a BGK scheme for recovering the pre collision PDF value.
         name: optional name of the boundary.
         init_wall_distance: Python callback function to calculate the wall distance for each cell near to the  boundary
         data_type: data type of the wall distance q
@@ -220,7 +221,7 @@ class QuadraticBounceBack(LbBoundary):
     def additional_data_init_callback(self):
         def default_callback(boundary_data, **_):
             for cell in boundary_data.index_array:
-                cell['q'] = -1
+                cell['q'] = 0.5
 
         if self.init_wall_distance:
             return self.init_wall_distance
@@ -238,6 +239,7 @@ class QuadraticBounceBack(LbBoundary):
         f_xf_inv = sp.Symbol("f_xf_inv")
         q = sp.Symbol("q")
         one = sp.Float(1.0)
+        two = sp.Float(2.0)
 
         subexpressions.append(Assignment(f_xf, f_out(dir_symbol)))
         subexpressions.append(Assignment(f_xf_inv, f_out(inv_dir[dir_symbol])))
@@ -253,17 +255,28 @@ class QuadraticBounceBack(LbBoundary):
 
         eq = lb_method.get_equilibrium_terms()
         cond = list()
+        # For the equilibrium values we calculate feq[Dir] + feq[inv_Dir] (equation E.4)
         for i in range(lb_method.stencil.Q):
-            subexpressions.append(Assignment(TypedSymbol(f"f_eq{i}", self.data_type), eq[i]))
-            if i > 0:
+            if i == 0:
+                subexpressions.append(Assignment(TypedSymbol("f_eq0", self.data_type), eq[0]))
+            else:
+                inv_idx = lb_method.stencil.index(lb_method.stencil.inverse_stencil_entries[i])
+                subexpressions.append(Assignment(TypedSymbol(f"f_eq{i}", self.data_type), eq[i] + eq[inv_idx]))
                 cond.append((TypedSymbol(f"f_eq{i}", self.data_type), sp.Eq(dir_symbol, i)))
+
         cond.append((TypedSymbol("f_eq0", self.data_type), True))
         subexpressions.append(Assignment(feq, sp.Piecewise(*cond)))
 
-        s1 = ((one - q) / (one + q)) * ((f_xf - feq) / (one - self.relaxation_rate) + feq)
-        s2 = ((q / (one + q)) * (f_xf + f_xf_inv))
+        # equation E.4 first summand
+        e41 = (f_xf_inv - f_xf) / two
+        # equation E.4 second summand
+        e42 = (f_xf_inv + f_xf - self.relaxation_rate * feq) / (two - two * self.relaxation_rate)
+        # equation E.3
+        fw = ((one - q) * (e41 + e42)) + q * f_xf_inv
+        # equation E.1
+        result = (one / (q + one)) * fw + (q / (q + one)) * f_xf
 
-        boundary_assignments = [Assignment(f_in(inv_dir[dir_symbol]), s1 + s2)]
+        boundary_assignments = [Assignment(f_in(inv_dir[dir_symbol]), result)]
         return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)