From 19155c925040c9efc3d3552adeccc875df0574ca Mon Sep 17 00:00:00 2001 From: Hesam Salehipour Date: Thu, 15 Feb 2024 14:00:41 -0500 Subject: [PATCH] fixed a subtle bug in interpolated bounce back BC --- examples/CFD/cylinder2d.py | 12 ++++++++---- src/boundary_conditions.py | 9 ++++++--- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/examples/CFD/cylinder2d.py b/examples/CFD/cylinder2d.py index e7fab29..2c9887d 100644 --- a/examples/CFD/cylinder2d.py +++ b/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) diff --git a/src/boundary_conditions.py b/src/boundary_conditions.py index a2d0087..6d3352a 100644 --- a/src/boundary_conditions.py +++ b/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,))