fixed a subtle bug in interpolated bounce back BC

examples/CFD/cylinder2d.py

@@ -16,6 +16,8 @@
 
 5. Visualization: The simulation outputs data in VTK format for visualization. It also generates images of the velocity field. The data can be visualized using software like ParaView.
 
+# To run type:
+nohup python3 examples/CFD/cylinder2d.py > logfile.log &
 """
 import os
 import json
@@ -79,7 +81,7 @@ def output_data(self, **kwargs):
         if timestep == 0:
             self.CL_max = 0.0
             self.CD_max = 0.0
-        if timestep > 0.8 * t_max:
+        if timestep > 0.5 * niter_max:
             # compute lift and drag over the cyliner
             cylinder = self.BCs[0]
             boundary_force = cylinder.momentum_exchange_force(kwargs['f_poststreaming'], kwargs['f_postcollision'])
@@ -95,7 +97,7 @@ def output_data(self, **kwargs):
             self.CL_max = max(self.CL_max, cl)
             self.CD_max = max(self.CD_max, cd)
             print('error= {:07.6f}, CL = {:07.6f}, CD = {:07.6f}'.format(err, cl, cd))
-            save_image(timestep, u)
+            # save_image(timestep, u)
 
 # Helper function to specify a parabolic poiseuille profile
 poiseuille_profile  = lambda x,x0,d,umax: np.maximum(0.,4.*umax/(d**2)*((x-x0)*d-(x-x0)**2))
@@ -132,9 +134,11 @@ def output_data(self, **kwargs):
             'print_info_rate': int(10000 / scale_factor),
             'return_fpost': True    # Need to retain fpost-collision for computation of lift and drag
         }
+        # characteristic time
+        tc = prescribed_vel/diam
+        niter_max = int(100//tc)
         sim = Cylinder(**kwargs)
-        t_max = int(1000000 / scale_factor)
-        sim.run(t_max)
+        sim.run(niter_max)
         CL_list.append(sim.CL_max)
         CD_list.append(sim.CD_max)
 

src/boundary_conditions.py

@@ -1058,16 +1058,19 @@ def set_proximity_ratio(self):
         -------
         None. The function updates the object's weights attribute in place.
         """
+        epsilon = 1e-12
+        nbd = len(self.indices[0])
         idx = np.array(self.indices).T
-        self.weights = np.full((idx.shape[0], self.lattice.q), 0.5)
+        bindex = np.arange(nbd)[:, None]
+        weights = np.full((idx.shape[0], self.lattice.q), 0.5)
         c = np.array(self.lattice.c)
         sdf_f = self.implicit_distances[self.indices]
         for q in range(1, self.lattice.q):
             solid_indices = idx + c[:, q]
             solid_indices_tuple = tuple(map(tuple, solid_indices.T))
             sdf_s = self.implicit_distances[solid_indices_tuple]
-            mask = self.iknownMask[:, q]
-            self.weights[mask, q] = sdf_f[mask] / (sdf_f[mask] - sdf_s[mask])
+            weights[:, q] = sdf_f / (sdf_f - sdf_s + epsilon)
+        self.weights = weights[bindex, self.iknown]
         return
 
     @partial(jit, static_argnums=(0,))