Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fixed a subtle bug in interpolated bounce back BC #36

Merged
merged 1 commit into from
Feb 20, 2024
Merged

fixed a subtle bug in interpolated bounce back BC #36

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 8 additions & 4 deletions examples/CFD/cylinder2d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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'])
Expand All @@ -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))
Expand Down Expand Up @@ -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)

Expand Down
9 changes: 6 additions & 3 deletions src/boundary_conditions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,))
Expand Down
Loading