19th Step of the 12 Steps to Navier-Stokes 😑 (3D-Lid Driven Cavity)

     One fine morning (in abundant spare time ⏳, of course), yours truly decided to code the 3D incompressible Navier–Stokes 🍃 equations using the finite-difference method. This post has the results of this adventure 🏞️ (so-far). As is customary with all my CFD work using commercial and home-made CFD codes, this too is an unofficial continuation of the series by Dr. Lorena Barba.

The code 🖳 shared in this post is fully vectorized 😲. Steps 13 to 18 are available here. 13, 14, 15, 16, 17, 18. Validation of the 2D code is available here and here. 🤓

     NOTE: This method requires a GPU, if dear readers don't have a GPU then please stop being peasants... 🙉

     If for some reason, you plan to use these codes in your scholarly work, do cite this blog as:

     Fahad Butt (2026). S-IBM (https://fluiddynamicscomputer.blogspot.com/2026/08/19th-step-of-12-steps-to-navier-stokes.html), Blogger. Retrieved Month Date, Year

     Lid-Driven Cavity 🕳 provides a complex 💢 flow 🌬 with very simple boundary conditions. Literally, everyone else uses this case to validate the code they write. The lid-driven cavity case is solved at ~Re 1,000 without any turbulence models or wall functions. The cavity used in the simulation is 1 x 1 x 1 m. This code is a simple 3D extension of this code. Well, the code is double in length as the 3D version requires under-relaxation and has one more momentum equation. Still, it is less than 100 lines ❗ The code is available here.

#%% import libraries

import cupy as cp

import matplotlib.pyplot as plt

#%% define parameters

l_cr = 1 # characteristic length

h = l_cr / 100 # grid spacing

dt = 0.001 # time step

L = 1 # domain length

D = 1 # domain depth

W = 1 # domain width

Nx = round(L / h) # grid points in x-axis

Ny = round(D / h) # grid points in y-axis

Nz = round(W / h) # grid points in z-axis

nu = 1 / 1000 # kinematic viscosity

Uinf = 1 # free stream velocity / inlet velocity / lid velocity

cfl = dt * Uinf / h # cfl number

travel = 10 # times the disturbance travels entire length of computational domain

TT = travel * L / Uinf # total time

ns = int(TT / dt) # number of time steps

Re = round(l_cr * Uinf / nu) # Reynolds number

#%% intialization

u = cp.zeros((Nx, Ny, Nz)) # x-velocity

v = cp.zeros((Nx, Ny, Nz)) # y-velocity

w = cp.zeros((Nx, Ny, Nz)) # z-velocity

p = cp.zeros((Nx, Ny, Nz)) # pressure

X, Y, Z = cp.meshgrid(cp.linspace(0, L, Nx), cp.linspace(0, D, Ny), cp.linspace(0, W, Nz), indexing = 'ij') # spatial grid

#%% pre calculate for speed

P1 = 1 / (2 * h * dt)

P2 = 1 / (4 * h * h)

P3 = h**2

P4 = 1 / 6

P5 = (2 / Re) * dt / h**2

P6 = dt / h

P7 = 1 - (6 * P5)

P8 = 0.75 # under relaxation

P9 = 1 - P8

#%% solve 3D Navier-Stokes equations

for nt in range(ns):

    dudx = u[2:, 1:-1, 1:-1] - u[:-2, 1:-1, 1:-1]

    dvdy = v[1:-1, 2:, 1:-1] - v[1:-1, :-2, 1:-1]

    dwdz = w[1:-1, 1:-1, 2:] - w[1:-1, 1:-1, :-2]

    pn = p.copy()

    b = P1 * (dudx + dvdy + dwdz) - P2 * (dudx**2 + dvdy**2 + dwdz**2 + 2 * ((u[1:-1, 2:, 1:-1] - u[1:-1, :-2, 1:-1]) * (v[2:, 1:-1, 1:-1] - v[:-2, 1:-1, 1:-1]) + (u[1:-1, 1:-1, 2:] - u[1:-1, 1:-1, :-2]) * (w[2:, 1:-1, 1:-1] - w[:-2, 1:-1, 1:-1]) + (v[1:-1, 1:-1, 2:] - v[1:-1, 1:-1, :-2]) * (w[1:-1, 2:, 1:-1] - w[1:-1, :-2, 1:-1]))) # divergence

    p[1:-1, 1:-1, 1:-1] = P4 * (pn[2:, 1:-1, 1:-1] + pn[:-2, 1:-1, 1:-1] + pn[1:-1, 2:, 1:-1] + pn[1:-1, :-2, 1:-1] + pn[1:-1, 1:-1, 2:] + pn[1:-1, 1:-1, :-2] - P3 * b) # mass

    p[0, :, :] = p[1, :, :] # dp/dx = 0 at x = 0

    p[-1, :, :] = p[-2, :, :] # dp/dx = 0 at x = L

    p[:, 0, :] = p[:, 1, :] # dp/dy = 0 at y = 0

    p[:, -1, :] = p[:, -2, :] # dp/dy = 0 at y = D

    p[:, :, 0] = p[:, :, 1] # dp/dz = 0 at z = 0

    p[:, :, -1] = p[:, :, -2] # dp/dz = 0 at z = H

    p = P8 * p + P9 * pn

    un = u.copy()

    vn = v.copy()

    wn = w.copy()

    u[1:-1, 1:-1, 1:-1] = un[1:-1, 1:-1, 1:-1] * P7 - P6 * (un[1:-1, 1:-1, 1:-1] * (un[2:, 1:-1, 1:-1] - un[:-2, 1:-1, 1:-1]) + vn[1:-1, 1:-1, 1:-1] * (un[1:-1, 2:, 1:-1] - un[1:-1, :-2, 1:-1]) + wn[1:-1, 1:-1, 1:-1] * (un[1:-1, 1:-1, 2:] - un[1:-1, 1:-1, :-2]) + p[2:, 1:-1, 1:-1] - p[:-2, 1:-1, 1:-1]) + P5 * (un[2:, 1:-1, 1:-1] + un[:-2, 1:-1, 1:-1] + un[1:-1, 2:, 1:-1] + un[1:-1, :-2, 1:-1] + un[1:-1, 1:-1, 2:] + un[1:-1, 1:-1, :-2]) # x momentum

    u[0, :, :] = 0 # u = 0 at x = 0

    u[-1, :, :] = 0 # u = 0 at # x = L

    u[:, 0, :] = 0 # u = 0 at y = 0

    u[:, -1, :] = 0 # u = 0 at y = D

    u[:, :, 0] = 0 # u = 0 at z = 0

    u[:, :, -1] = Uinf # u = Uinf at z = W

    u = P8 * u + P9 * un

    v[1:-1, 1:-1, 1:-1] = vn[1:-1, 1:-1, 1:-1] * P7 - P6 * (un[1:-1, 1:-1, 1:-1] * (vn[2:, 1:-1, 1:-1] - vn[:-2, 1:-1, 1:-1]) + vn[1:-1, 1:-1, 1:-1] * (vn[1:-1, 2:, 1:-1] - vn[1:-1, :-2, 1:-1]) + wn[1:-1, 1:-1, 1:-1] * (vn[1:-1, 1:-1, 2:] - vn[1:-1, 1:-1, :-2]) + p[1:-1, 2:, 1:-1] - p[1:-1, :-2, 1:-1]) + P5 * (vn[2:, 1:-1, 1:-1] + vn[:-2, 1:-1, 1:-1] + vn[1:-1, 2:, 1:-1] + vn[1:-1, :-2, 1:-1] + vn[1:-1, 1:-1, 2:] + vn[1:-1, 1:-1, :-2]) # y momentum

    v[0, :, :] = 0 # v = 0 at x = 0

    v[-1, :, :] = 0 # v = 0 at # x = L

    v[:, 0, :] = 0 # v = 0 at y = 0

    v[:, -1, :] = 0 # v = 0 at y = D

    v[:, :, 0] = 0 # v = 0 at z = 0

    v[:, :, -1] = 0 # v = 0 at z = W

    v = P8 * v + P9 * vn

    w[1:-1, 1:-1, 1:-1] = wn[1:-1, 1:-1, 1:-1] * P7 - P6 * (un[1:-1, 1:-1, 1:-1] * (wn[2:, 1:-1, 1:-1] - wn[:-2, 1:-1, 1:-1]) + vn[1:-1, 1:-1, 1:-1] * (wn[1:-1, 2:, 1:-1] - wn[1:-1, :-2, 1:-1]) + wn[1:-1, 1:-1, 1:-1] * (wn[1:-1, 1:-1, 2:] - wn[1:-1, 1:-1, :-2]) + p[1:-1, 1:-1, 2:] - p[1:-1, 1:-1, :-2]) + P5 * (wn[2:, 1:-1, 1:-1] + wn[:-2, 1:-1, 1:-1] + wn[1:-1, 2:, 1:-1] + wn[1:-1, :-2, 1:-1] + wn[1:-1, 1:-1, 2:] + wn[1:-1, 1:-1, :-2]) # z momentum

    w[0, :, :] = 0 # w = 0 at x = 0

    w[-1, :, :] = 0 # w = 0 at # x = L

    w[:, 0, :] = 0 # w = 0 at y = 0

    w[:, -1, :] = 0 # w = 0 at y = D

    w[:, :, 0] = 0 # w = 0 at z = 0

    w[:, :, -1] = 0 # w = 0 at z = W

    w = P8 * w + P9 * wn

    if nt % 1000 == 0:

        print("% Complete", 100 * nt / ns)

#%% post process

fig = plt.figure(dpi = 500)

ax = fig.add_subplot(111, projection = '3d')

ax.contourf(X[:, Ny // 2, :].get(), u[:, Ny // 2, :].get(), Z[:, Ny // 2, :].get(), zdir = 'y', offset = D / 2, levels = 128, cmap='jet', alpha = 0.5) # vertical plane

ax.set_xlim(0, L)

ax.set_ylim(0, D)

ax.set_zlim(0, W)

ax.set_xticks([0, L])

ax.set_yticks([0, D])

ax.set_zticks([0, W])

ax.tick_params(axis='x', pad = -2)

ax.tick_params(axis='y', pad = -2)

ax.tick_params(axis='z', pad = -2)

ax.set_xlabel('x [m]', labelpad = -15)

ax.set_ylabel('y [m]', labelpad = -15)

ax.set_zlabel('z [m]', labelpad = -15)

plt.gca().set_aspect('equal')

ax.view_init(elev = 22.5, azim = -45)

plt.show()

     The results from post processing are shown in Fig. 1. Within Fig. 1, the u, v and w components of velocities are shown along the plane of flow. The resulting pressure field and the velocity streamlines are shown within Fig. 2.

Fig. 1, Velocity components

Fig. 2, The pressure field and streamlines

     If you want to hire me as your next shining post-doc or collaborate in research, please reach out! Thank you very much for reading!

A Simple Poisson's Equation for Pressure

     In abundant spare time 🕰️, yours truly has removed ❌ the non-linear 〰️ terms from the Pressure Poisson Equation (PPE) which is derived 🎡 by summing the divergence of momentum equations and then applying conservation of mass ⚖️. The non-linear terms create a problem in convergence. Therefore, inspired by the style of work of the "Skipper"🐧, yours truly removed the problem causing non-linear terms 😀. The derived PPE is mentioned by equation 1. The PPE used in the code yours truly has developed is mention in equation 2. Within equations 1 and 2, the u and v are components of velocity along x and y-axis. The pressure is represented by p and while, t represents the time.

∂²p/∂x² + ∂²p/∂y² = 1/∆t * (∂u/∂x + ∂v/∂y - (∂u/∂x)² - (∂v/∂y)² - 2*∂v/∂y*∂u/∂x) [1]

∂²p/∂x² + ∂²p/∂y² = 1/(2 * ∆t) * (∂u/∂x + ∂v/∂y) [2]

     For the same grid and for solving the same problem, the time-step supported by equation 2 is 40x more as compared to equation 1. Therefore, the resultant compute per watt is 40x less for equation 2 as compared to equation 1 🤯. The code for equation 1 is shown first, followed by the code for equation 2. The results are compared via streamlines and pressure contours, with in Fig. 1. The benchmark case solved is of the lid-driven cavity which offer simple implementation and very complex flow physics. For validation of the code, refer to here, here and here.

Copyright <2025> <FAHAD BUTT>

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

     If you are mad 👨‍🔬 enough to use this code in your scholarly 👨‍🏫 work, then please remember to cite: Fahad Butt (2025). S-PPE (https://fluiddynamicscomputer.blogspot.com/2025/11/a-simple-poissons-equation-for-pressure.html), Blogger. Retrieved Month Date, Year

     MANDATORY NOTE: This method requires a GPU, if dear readers don't have a GPU then please stop being peasants... 🙉

Code 01

#%% import libraries

import cupy as cp

import matplotlib.pyplot as plt

#%% define parameters

l_cr = 1 # characteristic length

h = l_cr / 100 # grid spacing

dt = 0.00005 # time step

L = 1 # domain length

D = 1 # domain depth

Nx = round(L / h) + 1 # grid points in x-axis

Ny = round(D / h) + 1 # grid points in y-axis

nu = 1 / 400 # kinematic viscosity

Uinf = 1 # free stream velocity / inlet velocity / lid velocity

cfl = dt * Uinf / h # cfl number

travel = 10 # times the disturbance travels entire length of computational domain

TT = travel * L / Uinf # total time

ns = int(TT / dt) # number of time steps

Re = round(l_cr * Uinf / nu) # Osborne Reynolds and his number :)

#%% intialization (t = 0)

u = cp.zeros((Nx, Ny)) # x-velocity

v = cp.zeros((Nx, Ny)) # y-velocity

p = cp.zeros((Nx, Ny)) # pressure

X, Y = cp.meshgrid(cp.linspace(0, L, Nx), cp.linspace(0, D, Ny), indexing = 'ij') # spatial grid

#%% pre calculate for speed

P1 = h / (8 * dt)

P2 = (2 / Re) * dt / h**2

P3 = dt / h

P4 = 1 - (4 * P2)

#%% solve 2D incompressible Navier-Stokes equations

for _ in range(ns):

    pn = p.copy()

    p[1:-1, 1:-1] = 0.25 * (pn[2:, 1:-1] + pn[:-2, 1:-1] + pn[1:-1, 2:] + pn[1:-1, :-2]) - P1 * ((u[2:, 1:-1] - u[:-2, 1:-1] + v[1:-1, 2:] - v[1:-1, :-2]) - (u[2:, 1:-1] - u[:-2, 1:-1])**2 - (v[1:-1, 2:] - v[1:-1, :-2])**2 - (2 * (u[2:, 1:-1] - u[:-2, 1:-1]) * (v[1:-1, 2:] - v[1:-1, :-2]))) # pressure

    p[0, :] = p[1, :] # dp/dx = 0 at x = 0

    p[-1, :] = p[-2, :] # dp/dx = 0 at x = L

    p[:, 0] = p[:, 1] # dp/dy = 0 at y = 0

    p[:, -1] = p[:, -2] # dp/dy = 0 at y = D

    un = u.copy()

    vn = v.copy()

    u[1:-1, 1:-1] = un[1:-1, 1:-1] * P4 - P3 * (un[1:-1, 1:-1] * (un[2:, 1:-1] - un[:-2, 1:-1]) + vn[1:-1, 1:-1] * (un[1:-1, 2:] - un[1:-1, :-2]) + p[2:, 1:-1] - p[:-2, 1:-1]) + P2 * (un[2:, 1:-1] + un[:-2, 1:-1] + un[1:-1, 2:] + un[1:-1, :-2]) # x momentum

    u[0, :] = 0 # u = Uinf at x = 0

    u[-1, :] = 0 # u = 0 at x = L

    u[:, 0] = 0 # u = 0 at y = 0

    u[:, -1] = Uinf # u = Uinf at y = D

    v[1:-1, 1:-1] = vn[1:-1, 1:-1] * P4 - P3 * (un[1:-1, 1:-1] * (vn[2:, 1:-1] - vn[:-2, 1:-1]) + vn[1:-1, 1:-1] * (vn[1:-1, 2:] - vn[1:-1, :-2]) + p[1:-1, 2:] - p[1:-1, :-2]) + P2 * (vn[2:, 1:-1] + vn[:-2, 1:-1] + vn[1:-1, 2:] + vn[1:-1, :-2]) # y momentum

    v[0, :] = 0 # v = 0 at x = 0

    v[-1, :] = 0 # v = 0 at x = L

    v[:, 0] = 0 # v = 0 at y = 0

    v[:, -1] = 0 # v = 0 at y = D

Code 2

#%% import libraries

import cupy as cp

import matplotlib.pyplot as plt

#%% define parameters

l_cr = 1 # characteristic length

h = l_cr / 100 # grid spacing

dt = 0.002 # time step

L = 1 # domain length

D = 1 # domain depth

Nx = round(L / h) + 1 # grid points in x-axis

Ny = round(D / h) + 1 # grid points in y-axis

nu = 1 / 400 # kinematic viscosity

Uinf = 1 # free stream velocity / inlet velocity / lid velocity

cfl = dt * Uinf / h # cfl number

travel = 10 # times the disturbance travels entire length of computational domain

TT = travel * L / Uinf # total time

ns = int(TT / dt) # number of time steps

Re = round(l_cr * Uinf / nu) # Osborne Reynolds and his number :)

#%% intialization (t = 0)

u = cp.zeros((Nx, Ny)) # x-velocity

v = cp.zeros((Nx, Ny)) # y-velocity

p = cp.zeros((Nx, Ny)) # pressure

X, Y = cp.meshgrid(cp.linspace(0, L, Nx), cp.linspace(0, D, Ny), indexing = 'ij') # spatial grid

#%% pre calculate for speed

P1 = h / (16 * dt)

P2 = (2 / Re) * dt / h**2

P3 = dt / h

P4 = 1 - (4 * P2)

#%% solve 2D incompressible Navier-Stokes equations

for _ in range(ns):

    pn = p.copy()

    p[1:-1, 1:-1] = 0.25 * (pn[2:, 1:-1] + pn[:-2, 1:-1] + pn[1:-1, 2:] + pn[1:-1, :-2]) - P1 * (u[2:, 1:-1] - u[:-2, 1:-1] + v[1:-1, 2:] - v[1:-1, :-2]) # mass

    p[0, :] = p[1, :] # dp/dx = 0 at x = 0

    p[-1, :] = p[-2, :] # dp/dx = 0 at x = L

    p[:, 0] = p[:, 1] # dp/dy = 0 at y = 0

    p[:, -1] = p[:, -2] # dp/dy = 0 at y = D

    un = u.copy()

    vn = v.copy()

    u[1:-1, 1:-1] = un[1:-1, 1:-1] * P4 - P3 * (un[1:-1, 1:-1] * (un[2:, 1:-1] - un[:-2, 1:-1]) + vn[1:-1, 1:-1] * (un[1:-1, 2:] - un[1:-1, :-2]) + p[2:, 1:-1] - p[:-2, 1:-1]) + P2 * (un[2:, 1:-1] + un[:-2, 1:-1] + un[1:-1, 2:] + un[1:-1, :-2]) # x momentum

    u[0, :] = 0 # u = Uinf at x = 0

    u[-1, :] = 0 # u = 0 at x = L

    u[:, 0] = 0 # u = 0 at y = 0

    u[:, -1] = Uinf # u = Uinf at y = D

    v[1:-1, 1:-1] = vn[1:-1, 1:-1] * P4 - P3 * (un[1:-1, 1:-1] * (vn[2:, 1:-1] - vn[:-2, 1:-1]) + vn[1:-1, 1:-1] * (vn[1:-1, 2:] - vn[1:-1, :-2]) + p[1:-1, 2:] - p[1:-1, :-2]) + P2 * (vn[2:, 1:-1] + vn[:-2, 1:-1] + vn[1:-1, 2:] + vn[1:-1, :-2]) # y momentum

    v[0, :] = 0 # v = 0 at x = 0

    v[-1, :] = 0 # v = 0 at x = L

    v[:, 0] = 0 # v = 0 at y = 0

    v[:, -1] = 0 # v = 0 at y = D

Fig. 1, pressure range for both cases is 0 (blue) till 1 (red).

          If you want to hire me as your next shining post-doc or collaborate in research, please reach out! Thank you for reading!

Saithe Fish Simulation: ANSYS Fluent Dynamic Mesh Setup

     One of the most famous post on the blog can be read here. Worryingly😁, many fellow researchers and readers are interested in the aerodynamics of flexible robots 🤓. In this post, the dynamic mesh 🕸 settings used are shared 🥰. These settings are used to reproduce 🖨️ the results from [1], all those years ago. All in a hope that this post helps the readers in their scholarly work! 🎩

     Once the UDF 💻 has been acquired, the next step is to apply the UDF to the airfoil 🐠 geometry correctly ✔️. The airfoil geometry at the first time-step 🕰 i.e. at t = 0 for UDF 02 obtained from [1] is made available here. Once on the dynamic 🏎️ mesh page, select the options shown in Fig. 1. The options selected in Fig. 1 show the default parameters. Within Fig. 1, "wing" refers to the named selection that includes the only the airfoil geometry. Named selections can be created during the meshing process. The "wing" named selection is shown in Fig. 3.

Fig. 1, The dynamic mesh settings

     Before following the settings in Fig. 1, do remember to compile the UDF. To compile the UDF, please use the settings shown in Fig. 2. After selecting the UDF, select the options as shown in the Fig. 2 and then select Build and Load.

Fig. 2, Compile UDF

Fig. 3, Named selection for the dynamic mesh

     The maximum Lift ⬆️ force coefficient from the simulations performed using the method explained here is at 1.77 as compared to 1.68 [1]. The average Drag ⬅️ coefficient is at 0.097 as compared to 0.103 [1]. The obtained flow-field 🏟️ is shown in Fig. 4. Within Fig. 4, top row has v and u components of velocity while the bottom row shows pressure field❗

Fig. 4, The flow-field

     If you are still having trouble 😟, switch to immersed 🛀 boundary method. The immersed boundary method code yours truly wrote 🤓, is available here. Of course, this was done in abundant spare time 🕰️. The validation of this code is available here, here, here and more generally here 😼.

     If you want to hire me as your next shining post-doc or collaborate in research, please reach out! Thank you for reading!

References

[1] Shi, Fulong, Xin, Jianjian and Ou, Chuanzhong, Li, Zhiwei, Chang, Xing, Wan, Ling, "Effects of the Reynolds number and attack angle on wake dynamics of fish swimming in oblique flows", Physics of Fluids, 37(2), 025205, 2025 https://doi.org/10.1063/5.0252506 

A GPU Accelerated Simplified Immersed Boundary Method using Ray Casting

     Yours truly is an avid gamer, @fadoobaba (YouTube). The ray casting algorithm🛸; a fundamental 🧱 technique used in video game 🎮 development and computer graphics, has been implemented within the finite difference method 🏁 code yours truly has been developing. In abundant spare time, of course😼. In this post, this method is explained. For details about ray casting using matplotlib.path, refer to here. For validation of the code, refer to here, backwards-facing step, curved boundaries, here (moving cylinder) and here.

NOTE: This method requires a GPU, if dear readers don't have a GPU then please stop being peasants... 🙉

     If you plan to use these codes in your scholarly work, do cite this blog as:

     Fahad Butt (2025). S-IBM (https://fluiddynamicscomputer.blogspot.com/2025/07/a-simplified-immersed-boundary-method.html), Blogger. Retrieved Month Date, Year

     The first step is to setup the polygon 🐠. The polygon = cp.column_stack((x1, y1)) statement is used to combine x1 and y1 into a single array representing the polygon vertices. The following statements are used to store the x and y coordinates of the polygon vertices and nvert is the total number of vertices.

px = polygon[:, 0]

py = polygon[:, 1]

nvert = len(polygon)

     The boolean masks 👺 are then initialized. The following arrays track whether a grid point is inside the polygon based on horizontal and vertical ray intersections. False = outside ❌, True = inside ✅. Then a loop is implemented to for each edge 📏 of the polygon, defined by vertices i (current) and j (previous), with % nvert ensuring the polygon is closed (last vertex connects back to the first). (xi, yi) and (xj, yj) define the current edge and previous edge.

horizontal_inside = cp.zeros_like(test_x, dtype=bool) 

vertical_inside = cp.zeros_like(test_x, dtype=bool) 

for i in range(nvert):

    j = (i - 1) % nvert

xi, yi = px[i], py[i]

xj, yj = px[j], py[j]

     A check ⁉️ is performed to ensure a ray crosses the polygon edge once. cond1 statement checks if the test point lies between the y-values of the edge's endpoints i.e., a ray could cross it. intersect_x finds where the edge crosses a horizontal line at test_y. cond2 checks if the intersection lies to the right of the test point. ^= is XOR toggles the "inside" state each time the ray crosses an edge. The addition of small term ~1e-16 prevents division by zero for horizontal edges. Similar method is applied to verify the points using vertical ray. cond3 checks if the edge crosses a vertical ray (top to bottom) from (test_x, test_y). slope1 and intersect_y computes where the vertical ray at x=test_x intersects the edge. A point is inside only if both horizontal and vertical rays classify it as inside. A point is considered inside the foil only if both ray checks are true. curve is the 2D boolean mask of grid points inside the foil body.

cond1 = ((yi > test_y) != (yj > test_y))

slope = (xj - xi) / (yj - yi + 1e-16)

intersect_x = slope * (test_y - yi) + xi

cond2 = test_x < intersect_x 

horizontal_inside ^= cond1 & cond2

cond3 = ((xi > test_x) != (xj > test_x))

slope1 = (yj - yi) / (xj - xi + 1e-16)

intersect_y = slope1 * (test_x - xi) + yi

cond4 = test_y < intersect_y

vertical_inside ^= cond3 & cond4

inside = horizontal_inside & vertical_inside

curve = inside.reshape(X.T.shape)

     Using both horizontal and vertical rays reduces false positives (e.g., near sharp corners). The & ensures only points unambiguously inside are marked. Within Fig. 1, the points on the shape boundary, points inside 🪰 and outside 🐌 the shape boundary are shown. A point is a boundary if it is inside the body, but any of its neighbors is outside. Following statements are used to mark the boundary of the object. Boundary detects the "skin" of the body for applying no-slip, force, or stress conditions.

interior = curve[1:-1, 1:-1]

right = curve[2:, 1:-1]

left = curve[:-2, 1:-1]

top = curve[1:-1, 2:]

bottom = curve[1:-1, :-2]

boundary = interior & ~(right & left & top & bottom)

     The statements boundary_indices = cp.where(boundary_mask) ... valid_right = (boundary_indices[0] + 1 < X.shape[0]) & (~curve[right_neighbors]) ... etc. extract boundary indices and their valid neighbors i.e. these lines extract the (i, j) grid indices of the surface and locate which neighbor cells are valid fluid neighbors (outside the body and within domain bounds). These are needed to compute normals or apply boundary conditions via interpolation or extrapolation from the fluid.

Fig. 1, Mesh cells

     For validation of the results from present simulations, the case of flow around a circular cylinder is selected. Fig. 1 shows the results from the code at Re 200. The drag coefficient obtained from this code is 1.396 while from the literature, the value is at 1.4 [1]. The lift coefficient from the code is at 0.000134. Within Fig. 2, top row shows u and v velocity components and bottom row shows pressure and vorticity.

Fig. 2, The flow-field

     For the second an more complex validation case, a swimming fish is simulated. The lift and drag coefficients obtained from the simulation are compared with experimental results. The drag coefficient from the current code is at 0.342 while from the published literature, the value is at 0.348. Maximum lift coefficient is at 7.243 and 8 from the present code as compared to the published literature [2]. Within Fig. 3 the u and v velocities, pressure and velocity streamlines are shown for St = 0.8 and Reynolds number of 500. The computational mesh near the fish is shown in Fig. 4. Within Fig.4, a zoomed in view towards the right shows the mesh at the trailing edge of the fish.

Fig.3, Post-processing of results

Fig. 4, The mesh

     This method allows handling arbitrary deforming / non-deforming shapes on a fixed Cartesian grid. In summary, the method has the following steps.

1. Generate polygon shape (e.g., airfoil, cylinder)

2. Flatten mesh for vectorized testing

3. Use ray casting to check if points are inside shape

4. Build a boolean mask of body region

5. Identify surface (boundary) points

6. Extract boundary indices for physics coupling

References

     [1] Braza M, Chassaing P, Minh HH. Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder. Journal of Fluid Mechanics. 1986;165:79-130. doi:10.1017/S0022112086003014

     [2] Fulong Shi, Jianjian Xin, Chuanzhong Ou, Zhiwei Li, Xing Chang, Ling Wan; Effects of the Reynolds number and attack angle on wake dynamics of fish swimming in oblique flows. Physics of Fluids 1 February 2025; 37 (2): 025205 doi.org/10.1063/5.0252506

     Thank you for reading! If you want to hire me as a post-doc researcher in the fields of thermo-fluids and / or fracture mechanics, do reach out!