diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 2ec7d230..c52915af 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -20,6 +20,18 @@ jobs:
         uses: actions/checkout@v3
         with:
           fetch-depth: 1
+          
+      - name: Set up Python
+        uses: actions/setup-python@v4
+        with:
+          python-version: '3.x'
+
+      - name: Install Dependencies
+        run: |
+          pip install black
+
+      - name: Black Format Check
+        run: black --check .
 
       - name: Run Tests
         run: |
diff --git a/examples/smoke.py b/examples/smoke.py
index ede48a72..c84c4520 100644
--- a/examples/smoke.py
+++ b/examples/smoke.py
@@ -7,18 +7,20 @@
 
 import genesis as gs
 
+
 @ti.data_oriented
 class Jet(object):
     def __init__(
-            self,
-            world_center,
-            jet_radius,
-            orbit_radius,
-            orbit_radius_vel,
-            orbit_init_degree,
-            orbit_tau,
-            sub_orbit_radius,
-            sub_orbit_tau):
+        self,
+        world_center,
+        jet_radius,
+        orbit_radius,
+        orbit_radius_vel,
+        orbit_init_degree,
+        orbit_tau,
+        sub_orbit_radius,
+        sub_orbit_tau,
+    ):
         self.world_center = ti.Vector(world_center)
         self.orbit_radius = orbit_radius
         self.orbit_radius_vel = orbit_radius_vel
@@ -35,8 +37,7 @@ def __init__(
     @ti.func
     def get_pos(self, t: float):
         rel_pos = ti.Vector([self.orbit_radius + t * self.orbit_radius_vel, 0.0, 0.0])
-        rot_mat = ti.math.rot_by_axis(
-            ti.Vector([0.0, 1.0, 0.0]), self.orbit_init_radian + t * self.orbit_tau)[:3, :3]
+        rot_mat = ti.math.rot_by_axis(ti.Vector([0.0, 1.0, 0.0]), self.orbit_init_radian + t * self.orbit_tau)[:3, :3]
         rel_pos = rot_mat @ rel_pos
         return rel_pos
 
@@ -73,9 +74,9 @@ def get_tan_dir(self, t: float):
 def main():
 
     ########################## init ##########################
-    gs.init(seed=0, precision='32', logging_level='debug')
+    gs.init(seed=0, precision="32", logging_level="debug")
 
-    video_path = Path(__file__).parent / 'video'
+    video_path = Path(__file__).parent / "video"
     video_path.mkdir(exist_ok=True, parents=True)
 
     res = 384
@@ -109,7 +110,9 @@ def main():
             orbit_tau=orbit_tau,
             sub_orbit_radius=sub_orbit_radius,
             jet_radius=jet_radius,
-            sub_orbit_tau=sub_orbit_tau) for orbit_init_degree in np.linspace(0, 360, 3, endpoint=False)
+            sub_orbit_tau=sub_orbit_tau,
+        )
+        for orbit_init_degree in np.linspace(0, 360, 3, endpoint=False)
     ]
     scene.sim.solvers[-1].set_jets(jets)
 
@@ -122,16 +125,16 @@ def main():
 
     for i in range(num_steps):
 
-        scalars = scene.sim.solvers[-1].grid.q.to_numpy().astype(np.float32) # (res, res, res, 3)
+        scalars = scene.sim.solvers[-1].grid.q.to_numpy().astype(np.float32)  # (res, res, res, 3)
         scalars[scalars < 1e-4] = 0
         layer = scalars[:, res // 2, :]
 
         rgb = (255 * layer).astype(np.uint8)
-        cv2.imwrite(str(video_path / f'{i:04d}.png'), cv2.cvtColor(rgb, cv2.COLOR_RGB2BGR))
+        cv2.imwrite(str(video_path / f"{i:04d}.png"), cv2.cvtColor(rgb, cv2.COLOR_RGB2BGR))
 
         for _ in range(substeps):
             scene.step()
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()
diff --git a/genesis/engine/entities/sf_entity.py b/genesis/engine/entities/sf_entity.py
index b349699a..fb0c8474 100644
--- a/genesis/engine/entities/sf_entity.py
+++ b/genesis/engine/entities/sf_entity.py
@@ -12,11 +12,14 @@
 
 @ti.data_oriented
 class SFParticleEntity(ParticleEntity):
-    '''
+    """
     PBD-based entity represented solely by particles.
-    '''
+    """
+
     def __init__(self, scene, solver, material, morph, surface, particle_size, idx, particle_start):
-        super().__init__(scene, solver, material, morph, surface, particle_size, idx, particle_start, need_skinning=False)
+        super().__init__(
+            scene, solver, material, morph, surface, particle_size, idx, particle_start, need_skinning=False
+        )
 
     def process_input(self, in_backward=False):
         # TODO: implement this
@@ -31,5 +34,5 @@ def sample(self):
         pass
 
     def update_particles(self, particles):
-        self._particles   = particles
+        self._particles = particles
         self._n_particles = particles.shape[0]
diff --git a/genesis/engine/materials/SF/__init__.py b/genesis/engine/materials/SF/__init__.py
index 2396c2e1..218f8cd8 100644
--- a/genesis/engine/materials/SF/__init__.py
+++ b/genesis/engine/materials/SF/__init__.py
@@ -1,2 +1,2 @@
 from .base import Base
-from .smoke import Smoke
\ No newline at end of file
+from .smoke import Smoke
diff --git a/genesis/engine/materials/SF/smoke.py b/genesis/engine/materials/SF/smoke.py
index 6210db0e..ef9f519f 100644
--- a/genesis/engine/materials/SF/smoke.py
+++ b/genesis/engine/materials/SF/smoke.py
@@ -10,4 +10,4 @@ def __init__(self):
 
     @property
     def sampler(self):
-        return 'regular'
+        return "regular"
diff --git a/genesis/engine/scene.py b/genesis/engine/scene.py
index 7856673f..8d4a248d 100644
--- a/genesis/engine/scene.py
+++ b/genesis/engine/scene.py
@@ -306,10 +306,12 @@ def add_entity(
 
         elif isinstance(material, (gs.materials.SF.Smoke)):
             if surface.vis_mode is None:
-                surface.vis_mode = 'particle'
+                surface.vis_mode = "particle"
 
-            if surface.vis_mode not in ['particle']:
-                gs.raise_exception(f"Unsupported `surface.vis_mode` for material {material}: '{surface.vis_mode}'. Expected one of: ['particle', 'recon'].")
+            if surface.vis_mode not in ["particle"]:
+                gs.raise_exception(
+                    f"Unsupported `surface.vis_mode` for material {material}: '{surface.vis_mode}'. Expected one of: ['particle', 'recon']."
+                )
 
         elif isinstance(material, (gs.materials.PBD.Base, gs.materials.MPM.Base, gs.materials.SPH.Base)):
             if surface.vis_mode is None:
diff --git a/genesis/engine/solvers/sf_solver.py b/genesis/engine/solvers/sf_solver.py
index 79ffdf73..e271876b 100644
--- a/genesis/engine/solvers/sf_solver.py
+++ b/genesis/engine/solvers/sf_solver.py
@@ -12,9 +12,9 @@
 
 @ti.data_oriented
 class SFSolver(Solver):
-    '''
+    """
     Stable Fluid solver for eulerian-based gaseous simulation.
-    '''
+    """
 
     def __init__(self, scene, sim, options):
         super().__init__(scene, sim, options)
@@ -22,14 +22,14 @@ def __init__(self, scene, sim, options):
         if options is None:
             return
 
-        self.n_grid        = options.res
-        self.dx            = 1 / self.n_grid
-        self.res           = (self.n_grid, self.n_grid, self.n_grid)
-        self.solver_iters  = options.solver_iters
-        self.decay         = options.decay
+        self.n_grid = options.res
+        self.dx = 1 / self.n_grid
+        self.res = (self.n_grid, self.n_grid, self.n_grid)
+        self.solver_iters = options.solver_iters
+        self.decay = options.decay
 
-        self.t             = 0.0
-        self.inlet_s       = options.inlet_s
+        self.t = 0.0
+        self.inlet_s = options.inlet_s
 
         self.jets = []
 
@@ -48,22 +48,21 @@ def is_active(self):
 
     def setup_fields(self):
         cell_state = ti.types.struct(
-            v     = gs.ti_vec3,
-            v_tmp = gs.ti_vec3,
-            div   = gs.ti_float,
-            p     = gs.ti_float,
-            q     = ti.types.vector(len(self.jets), gs.ti_float),
+            v=gs.ti_vec3,
+            v_tmp=gs.ti_vec3,
+            div=gs.ti_float,
+            p=gs.ti_float,
+            q=ti.types.vector(len(self.jets), gs.ti_float),
         )
 
-        self.grid    = cell_state.field(shape=self.res, layout=ti.Layout.SOA)
+        self.grid = cell_state.field(shape=self.res, layout=ti.Layout.SOA)
 
         # swap area for pressure projection solver
         self.p_swap = TexPair(
-            cur = ti.field(dtype=gs.ti_float, shape=self.res),
-            nxt = ti.field(dtype=gs.ti_float, shape=self.res),
+            cur=ti.field(dtype=gs.ti_float, shape=self.res),
+            nxt=ti.field(dtype=gs.ti_float, shape=self.res),
         )
 
-
     @ti.kernel
     def init_fields(self):
         for i, j, k in ti.ndrange(*self.res):
@@ -183,12 +182,11 @@ def compute_location(self, u, v, w, du, dv, dw):
         I = max(0, min(self.n_grid - 1, I))
         return I
 
-
     @ti.func
     def is_free(self, u, v, w, du, dv, dw):
         flag = 1
 
-        I = ti.Vector([int(u+du), int(v+dv), int(w+dw)])
+        I = ti.Vector([int(u + du), int(v + dv), int(w + dw)])
         if (I < 0).any() or (I > self.n_grid - 1).any():
             flag = 0
 
@@ -196,10 +194,10 @@ def is_free(self, u, v, w, du, dv, dw):
 
     @ti.func
     def trilerp_scalar(self, qf, p, qf_idx):
-        '''
+        """
         p: position, within (0, 1).
         qf: field for interpolation
-        '''
+        """
         # convert position to grid index
         base_I = ti.floor(p - 0.5, gs.ti_int)
         p_I = p - 0.5
@@ -219,10 +217,10 @@ def trilerp_scalar(self, qf, p, qf_idx):
 
     @ti.func
     def trilerp(self, qf, p):
-        '''
+        """
         p: position, within (0, 1).
         qf: field for interpolation
-        '''
+        """
         # convert position to grid index
         base_I = ti.floor(p - 0.5, gs.ti_int)
         p_I = p - 0.5
@@ -243,9 +241,9 @@ def trilerp(self, qf, p):
     # RK3
     @ti.func
     def backtrace(self, vf, p, dt):
-        '''
+        """
         vf: velocity field
-        '''
+        """
         v1 = self.trilerp(vf, p)
         p1 = p - 0.5 * dt * v1
         v2 = self.trilerp(vf, p1)
@@ -269,6 +267,7 @@ def set_state(self, f, state):
     def reset_grad(self):
         return None
 
+
 class TexPair:
     def __init__(self, cur, nxt):
         self.cur = cur
diff --git a/genesis/options/solvers.py b/genesis/options/solvers.py
index c82515de..54e80b18 100644
--- a/genesis/options/solvers.py
+++ b/genesis/options/solvers.py
@@ -483,15 +483,16 @@ class SFOptions(Options):
     dt : float, optional
         Time duration for each simulation step in seconds. If none, it will inherit from `SimOptions`. Defaults to None.
     """
-    dt           : Optional[float] = None
-    res          : Optional[int] = 128
-    solver_iters : Optional[int] = 500
-    decay        : Optional[float] = 0.99
-
-    T_low        : Optional[float] = 1.0
-    T_high       : Optional[float] = 0.0
-
-    inlet_pos    : Optional[tuple[int, int, int]] = (0.6, 0.0, 0.1)
-    inlet_vel    : Optional[tuple[int, int, int]] = (0, 0, 1)
-    inlet_quat   : Optional[tuple[int, int, int, int]] = (1, 0, 0, 0)
-    inlet_s      : Optional[float] = 400.0
+
+    dt: Optional[float] = None
+    res: Optional[int] = 128
+    solver_iters: Optional[int] = 500
+    decay: Optional[float] = 0.99
+
+    T_low: Optional[float] = 1.0
+    T_high: Optional[float] = 0.0
+
+    inlet_pos: Optional[tuple[int, int, int]] = (0.6, 0.0, 0.1)
+    inlet_vel: Optional[tuple[int, int, int]] = (0, 0, 1)
+    inlet_quat: Optional[tuple[int, int, int, int]] = (1, 0, 0, 0)
+    inlet_s: Optional[float] = 400.0
diff --git a/pyproject.toml b/pyproject.toml
index 4f0f3371..4dd212bb 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -34,6 +34,7 @@ dependencies = [
     "pymeshlab",
     "coacd",
     "OpenEXR",
+    "black",
 ]
 
 [tool.setuptools.packages.find]
@@ -48,6 +49,9 @@ genesis = [
     "ext/VolumeSampling",
 ]
 
+[tool.black]
+line-length = 120
+
 [project.scripts]
 gs = "genesis._main:main"