Chapter 8: Programming Practices#

Topics Covered:

  • Code documentation and commenting

  • Testing and validation

  • Error handling and exceptions

  • Debugging techniques

8.1 Code Documentation and Commenting#

Good documentation makes code understandable, maintainable, and reusable. In chemical engineering, where calculations can be complex and involve many parameters, clear documentation is critical.

8.1.1 Docstrings#

Docstrings are special comments that describe what a function does, its parameters, and return values.

import math

def reynolds_number(rho, v, D, mu):
    """
    Calculate the Reynolds number for fluid flow in a pipe.
    
    The Reynolds number is a dimensionless quantity that predicts
    flow patterns in different fluid flow situations.
    
    Parameters
    ----------
    rho : float
        Fluid density (kg/m³)
    v : float
        Flow velocity (m/s)
    D : float
        Pipe diameter (m)
    mu : float
        Dynamic viscosity (Pa·s)
    
    Returns
    -------
    float
        Reynolds number (dimensionless)
    
    Examples
    --------
    >>> reynolds_number(1000, 2.0, 0.1, 0.001)
    200000.0
    
    Notes
    -----
    Flow regimes:
    - Re < 2300: Laminar flow
    - 2300 ≤ Re ≤ 4000: Transitional flow
    - Re > 4000: Turbulent flow
    """
    Re = (rho * v * D) / mu
    return Re

# Test the function
Re = reynolds_number(1000, 2.0, 0.1, 0.001)
print(f"Reynolds number: {Re:.0f}")

# # Access the docstring
# print("\nFunction documentation:")
# print(reynolds_number.__doc__)
Reynolds number: 200000
print(reynolds_number.__doc__)
    Calculate the Reynolds number for fluid flow in a pipe.
    
    The Reynolds number is a dimensionless quantity that predicts
    flow patterns in different fluid flow situations.
    
    Parameters
    ----------
    rho : float
        Fluid density (kg/m³)
    v : float
        Flow velocity (m/s)
    D : float
        Pipe diameter (m)
    mu : float
        Dynamic viscosity (Pa·s)
    
    Returns
    -------
    float
        Reynolds number (dimensionless)
    
    Examples
    --------
    >>> reynolds_number(1000, 2.0, 0.1, 0.001)
    200000.0
    
    Notes
    -----
    Flow regimes:
    - Re < 2300: Laminar flow
    - 2300 ≤ Re ≤ 4000: Transitional flow
    - Re > 4000: Turbulent flow
    
import math, os


print(math.sqrt.__doc__)
print(math.sin.__doc__)

print(os.path.__doc__)
Return the square root of x.
Return the sine of x (measured in radians).
Common operations on Posix pathnames.

Instead of importing this module directly, import os and refer to
this module as os.path.  The "os.path" name is an alias for this
module on Posix systems; on other systems (e.g. Windows),
os.path provides the same operations in a manner specific to that
platform, and is an alias to another module (e.g. ntpath).

Some of this can actually be useful on non-Posix systems too, e.g.
for manipulation of the pathname component of URLs.

8.1.2 Inline Comments#

Use comments to explain why you’re doing something, not what you’re doing (the code shows what).

def arrhenius_rate_constant(A, Ea, T):
    """
    Calculate reaction rate constant using Arrhenius equation.
    
    Parameters
    ----------
    A : float
        Pre-exponential factor (1/s)
    Ea : float
        Activation energy (J/mol)
    T : float
        Temperature (K)
    
    Returns
    -------
    float
        Rate constant k (1/s)
    """
    R = 8.314  # J/(mol·K) - Universal gas constant
    
    # Negative sign because higher Ea means slower reaction
    exponent = -Ea / (R * T)
    
    k = A * math.exp(exponent)
    
    return k

# Good comment: explains WHY we use these specific values
# Using typical values for a first-order decomposition reaction
A = 1.0e10  # s⁻¹
Ea = 50000  # J/mol
T = 298.15  # K (25°C, standard conditions)

k = arrhenius_rate_constant(A, Ea, T)
print(f"Rate constant at {T} K: {k:.4e} s⁻¹")
Rate constant at 298.15 K: 1.7374e+01 s⁻¹

8.1.3 Code Organization Best Practices#

# BAD: Unclear variable names, no documentation
def calc(x, y, z):
    return x * y / z**2

# GOOD: Clear names, documentation, units specified
def gravitational_force(mass1_kg, mass2_kg, distance_m):
    """
    Calculate gravitational force between two masses.
    
    Parameters
    ----------
    mass1_kg : float
        Mass of first object (kg)
    mass2_kg : float
        Mass of second object (kg)
    distance_m : float
        Distance between centers (m)
    
    Returns
    -------
    float
        Gravitational force (N)
    """
    G = 6.674e-11  # N·m²/kg² - Gravitational constant
    
    force_N = G * mass1_kg * mass2_kg / (distance_m ** 2)
    
    return force_N

# Test both versions
F1 = calc(100, 50, 1.0)
F2 = gravitational_force(100, 50, 1.0)

print(f"Both give same result: {F1:.6e} N")
print("But the second is much more readable!")
gravitational_force.__doc__
Both give same result: 5.000000e+03 N
But the second is much more readable!
'\n    Calculate gravitational force between two masses.\n    \n    Parameters\n    ----------\n    mass1_kg : float\n        Mass of first object (kg)\n    mass2_kg : float\n        Mass of second object (kg)\n    distance_m : float\n        Distance between centers (m)\n    \n    Returns\n    -------\n    float\n        Gravitational force (N)\n    '

8.2 (Optional) Testing and Validation#

Testing ensures your code produces correct results. For chemical engineering calculations, errors can have serious consequences.

8.2.1 Unit Testing with Known Results#

def ideal_gas_pressure(n, T, V):
    """
    Calculate pressure using ideal gas law: P = nRT/V
    
    Parameters
    ----------
    n : float
        Amount of gas (mol)
    T : float
        Temperature (K)
    V : float
        Volume (m³)
    
    Returns
    -------
    float
        Pressure (Pa)
    """
    R = 8.314  # J/(mol·K)
    P = (n * R * T) / V
    return P

# Test with known result: 1 mol at STP should give ~101325 Pa
def test_ideal_gas_pressure():
    # Test 1: Known textbook example
    P1 = ideal_gas_pressure(n=50.0, T=373.15, V=0.5)
    print(f"Calculated Pressure: {P1:.2f} Pa")
    assert abs(P1 - 31097.41) < 1.0, "Test 1 failed"

    # Test 2: Room temperature gas
    P2 = ideal_gas_pressure(n=1.0, T=298.15, V=1.0)
    print(f"Calculated Pressure: {P2:.2f} Pa")
    assert abs(P2 - 2478.9) < 1.0, "Test 2 failed"

    # # Test 3: Scaling check (double moles → double pressure)
    # P3 = ideal_gas_pressure(n=2.0, T=300.0, V=1.0)
    # P4 = ideal_gas_pressure(n=1.0, T=300.0, V=1.0)
    # assert abs(P3 - 2 * P4) < 1e-6, "Test 3 failed"

    print("All ideal gas law tests passed!")


# Run tests
test_ideal_gas_pressure()
Calculated Pressure: 310236.91 Pa
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[6], line 44
     40     print("All ideal gas law tests passed!")
     43 # Run tests
---> 44 test_ideal_gas_pressure()

Cell In[6], line 28, in test_ideal_gas_pressure()
     26 P1 = ideal_gas_pressure(n=50.0, T=373.15, V=0.5)
     27 print(f"Calculated Pressure: {P1:.2f} Pa")
---> 28 assert abs(P1 - 31097.41) < 1.0, "Test 1 failed"
     30 # Test 2: Room temperature gas
     31 P2 = ideal_gas_pressure(n=1.0, T=298.15, V=1.0)

AssertionError: Test 1 failed

8.2.2 Boundary Testing#

def test_edge_cases():
    # Very small volume → large pressure
    P = ideal_gas_pressure(1.0, 300.0, 1e-3)
    assert P > 1e6

    # Zero moles → zero pressure
    P = ideal_gas_pressure(0.0, 300.0, 1.0)
    assert P == 0.0

    print("Edge case tests passed!")

# Run edge case tests
test_edge_cases()
Edge case tests passed!

(Optional) 8.3 Error Handling#

Robust code anticipates and handles errors gracefully. This prevents crashes and provides useful feedback.

8.3.1 Try-Except Blocks#

def safe_division(numerator, denominator):
    """
    Safely divide two numbers with error handling.
    
    Parameters
    ----------
    numerator : float
    denominator : float
    
    Returns
    -------
    float or None
        Result of division, or None if error
    """
    try:
        result = numerator / denominator
        return result
    except ZeroDivisionError:
        print(f"Error: Cannot divide {numerator} by zero!")
        return None
    except TypeError:
        print(f"Error: Invalid input types: {type(numerator)}, {type(denominator)}")
        return None

# Test error handling
print("Testing error handling:\n")

result1 = safe_division(10, 2)
print(f"10 / 2 = {result1}")

result2 = safe_division(10, 0)
print(f"10 / 0 = {result2}")

result3 = safe_division(10, "hello")
print(f"10 / 'hello' = {result3}")
Testing error handling:

10 / 2 = 5.0
Error: Cannot divide 10 by zero!
10 / 0 = None
Error: Invalid input types: <class 'int'>, <class 'str'>
10 / 'hello' = None

8.3.2 Input Validation#

def heat_capacity_calculation(mass, Cp, delta_T):
    """
    Calculate heat required: Q = m × Cp × ΔT
    
    Parameters
    ----------
    mass : float
        Mass (kg), must be positive
    Cp : float
        Specific heat capacity (J/kg·K), must be positive
    delta_T : float
        Temperature change (K), can be negative for cooling
    
    Returns
    -------
    float
        Heat energy (J)
    
    Raises
    ------
    ValueError
        If mass or Cp is not positive
    """
    # Validate inputs
    if mass <= 0:
        raise ValueError(f"Mass must be positive, got {mass} kg")
    
    if Cp <= 0:
        raise ValueError(f"Heat capacity must be positive, got {Cp} J/kg·K")
    
    # Calculate heat
    Q = mass * Cp * delta_T
    
    return Q
# Valid input
try:
    Q = heat_capacity_calculation(2.5, 4184, 50)
    print(f"✓ Valid: Q = {Q:.2f} J")
except ValueError as e:
    print(f"✗ {e}")
✓ Valid: Q = 523000.00 J
# Invalid mass
try:
    Q = heat_capacity_calculation(-2.5, 4184, 50)
    print(f"✓ Should not reach here")
except ValueError as e:
    print(f"✓ Caught error: {e}")
✓ Caught error: Mass must be positive, got -2.5 kg

8.3.3 Graceful Degradation#

def average_temperature(temperatures):
    """
    Calculate the average temperature from a list.
    Invalid entries are skipped.
    
    Parameters
    ----------
    temperatures : list
        List of temperature values (°C)
    
    Returns
    -------
    float
        Average temperature of valid values
    """
    total = 0
    count = 0

    for temp in temperatures:
        try:
            temp_value = float(temp)

            if temp_value < -273.15:
                raise ValueError("Below absolute zero")

            total += temp_value
            count += 1

        except (ValueError, TypeError):
            continue  # Skip invalid values

    if count == 0:
        return None

    return total / count
# Test data with valid and invalid values
temps = [22.5, 25.0, "hot", -300, 18.7, None, 20.1]

avg = average_temperature(temps)

if avg is not None:
    print(f"Average temperature: {avg:.2f} °C")
else:
    print("No valid temperature data found.")
Average temperature: 21.58 °C

8.4 Debugging Techniques#

Debugging is the process of finding and fixing errors in code. Systematic debugging saves time and frustration.

8.4.2 Assert Statements for Assumptions#

assert condition
def calculate_residence_time(volume, flow_rate):
    """
    Calculate residence time: τ = V / Q
    
    Parameters
    ----------
    volume : float
        Reactor volume (L)
    flow_rate : float
        Volumetric flow rate (L/min)
    
    Returns
    -------
    float
        Residence time (min)
    """
    # Assert our assumptions
    assert volume > 0, f"Volume must be positive, got {volume}"
    assert flow_rate > 0, f"Flow rate must be positive, got {flow_rate}"
    
    tau = volume / flow_rate
    
    # Check the result makes sense
    assert tau > 0, "Residence time should be positive"
    
    return tau

# Valid case
# try:
#     tau = calculate_residence_time(100, 5)
#     print(f"✓ Residence time: {tau} min")
# except AssertionError as e:
#     print(f"✗ Assertion failed: {e}")

# # Invalid case
# try:
#     tau = calculate_residence_time(-100, 5)
#     print(f"✓ Residence time: {tau} min")
# except AssertionError as e:
#     print(f"✓ Caught invalid input: {e}")
tau = calculate_residence_time(100, 5)
print(f"✓ Residence time: {tau} min")
✓ Residence time: 20.0 min

8.4.3 Usage of breakpoint()#

breakpoint() is a built-in Python function that pauses program execution at a specific line and enters interactive debugging mode. It allows you to inspect variable values, step through code line by line, and identify logic or runtime errors while the program is running. This is useful for understanding how a program behaves and for locating bugs without adding multiple print statements.

Run kelvin.py, and check how breakpoint() works.

Summary#

In this chapter, you learned:

  1. Documentation

    • Writing comprehensive docstrings

    • Using comments effectively

    • Naming conventions and code organization

  2. Testing

    • Unit testing with known results

    • Boundary condition testing

    • Conservation law verification

  3. Error Handling

    • Try-except blocks

    • Input validation

    • Graceful degradation

  4. Debugging

    • Print statement debugging

    • Assert statements

    • Step-by-step verification

  5. Best Practices

    • Professional code style

    • Comprehensive documentation

    • Robust error handling

    • Systematic testing

Remember: Good programming practices prevent errors, save debugging time, and make your code understandable to others (including your future self!).