Reading Time: 10 minutes

Key Takeaways

  • Debugging scientific Python code follows an escalation ladder: print, assert, pdb, profiling, and tests.
  • Most scientific bugs are numerical, such as NaN values, array shape mismatches, precision loss, or unstable time steps.
  • print() catches many beginner bugs, but assertions and debuggers make failures easier to isolate.
  • A reliable scientific debugging workflow is: make it fail reliably, isolate the problem, change one thing, verify, and add a test.
  • Profiling can reveal hidden bugs because slow code is often a sign of wrong array logic, accidental loops, or excessive allocation.

Debugging scientific Python code is one of the most frustrating skills to develop. You write what should be a simple heat diffusion simulation, the code runs without errors, and the output is completely wrong.

No crash. No traceback. Just silently incorrect numbers.

That is the unique challenge of scientific debugging. The most dangerous bugs are often the ones that do not crash at all. A wrong number can propagate through thousands of iterations and produce a result that looks physically plausible until you check conservation, stability, or convergence.

This guide explains how to debug scientific Python code at several levels: print statements, logging, assertions, interactive debugging, profiling, and regression tests.

Whether you are running a first FiPy simulation or debugging a production research pipeline, this workflow gives you a practical structure to follow.

What Is Different About Debugging Scientific Code?

Scientific debugging is different from debugging a web app, API, or ordinary data pipeline.

The numbers can lie. A web app often fails visibly: a button does not work, an API returns an error, or a database query fails. Scientific code can run perfectly and still produce wrong results. A negative concentration, exploding energy field, or hidden NaN may only appear if you check for it directly.

Correctness is defined by physics, not just syntax. A Python program that runs is syntactically valid. A scientific program is credible only if it matches known solutions, conserves expected quantities, respects boundary conditions, and converges under refinement.

Stack traces often point into NumPy, SciPy, or FiPy internals. The real bug may be in your setup, but the error message may appear deep inside a library. That is why scientific debugging needs a structured process, not just traceback reading.

The Debugging Workflow: A Five-Step Mental Model

A useful debugging workflow for scientific code follows five steps.

Step 1: Make It Fail Reliably

Your first job is to create a test case that fails every time. In scientific code, this usually means building a minimal reproducible example that produces the wrong result.

You do not need the full simulation. You need the smallest piece of code that demonstrates the bug.

import numpy as np

# Full simulation may be hundreds of lines.
# Start with the smallest failing update step.

phi = np.ones((10, 10))
phi[:] = 0.5

gradient = np.gradient(phi)[0]
phi_new = phi - 0.1 * gradient

print(f"min(phi_new) = {np.min(phi_new)}")
print(f"max(phi_new) = {np.max(phi_new)}")

By isolating the update step, you create a reproducible case. This is the foundation for every later debugging step.

Step 2: Divide and Conquer

Once you have a failing case, isolate it further. Find the module, function, line, or physics term responsible for the failure.

For simulation code, a useful tactic is to test each physics term separately. If your model combines diffusion, reaction, and advection, test each one alone before debugging the coupled system.

solution_diffusion = solve_diffusion_only(initial_conditions, dt)
assert np.all(solution_diffusion >= 0), "Diffusion produced negative values"

solution_reaction = solve_reaction_only(initial_conditions, dt)
assert np.all(solution_reaction >= 0), "Reaction produced negative values"

solution_coupled = solve_coupled(initial_conditions, dt)

This narrows the search space. Instead of debugging the whole simulation, you debug one term or one coupling mechanism.

Step 3: Change One Thing at a Time

Scientific debugging becomes confusing when you change several variables at once. If you change the time step, boundary conditions, solver tolerance, and initial condition together, you cannot know which change mattered.

Change one parameter. Run. Compare. Document.

Use version control to track changes. Even simple commits such as “Changed dt from 0.01 to 0.001” can save hours later.

Step 4: Use the Debugger

When prints and assertions are not enough, use pdb, ipdb, or an IDE debugger. A debugger lets you pause execution, inspect variables, step through code, and check why a value changed.

Step 5: Add to the Test Suite

After you fix a bug, add a test that would have caught it. This prevents the same bug from returning later.

def test_no_negative_phases():
    """Verify phase field stays non-negative after diffusion."""
    phi = np.ones((10, 10))
    phi = diffusion_step(phi, dt=0.01)
    assert np.min(phi) >= 0, f"Got min value {np.min(phi)}"

The workflow is simple: isolate, divide, change one thing, debug, and test. The rest of this guide explains which tools to use at each level.

Level 1: Print Statements and Logging

print() is still useful. It is the simplest debugging tool and catches many early bugs. The key is to use it with structure.

The Problem With Bad Prints

A common beginner mistake is printing values without context.

for i in range(100):
    print(phi[i])

After a few iterations, you have many numbers with no clear meaning. You do not know which variable, time step, or condition produced each line.

The Right Way: Structured Prints

Label everything. Include time step, variable name, minimum, maximum, and any conserved quantity you care about.

import numpy as np

print(
    f"[t={t:06.4f}] "
    f"min(phi)={np.min(phi):.6f}, "
    f"max(phi)={np.max(phi):.6f}"
)

print(
    f"[t={t:06.4f}] "
    f"energy={np.sum(phi * phi):.6f}, "
    f"mass={np.sum(phi):.6f}"
)

This output is easy to scan:

[t=0.0010] min(phi)=0.000000, max(phi)=1.000000
[t=0.0010] energy=0.850000, mass=1.000000
[t=0.0020] min(phi)=0.000000, max(phi)=1.023456
[t=0.0020] energy=0.873456, mass=1.012345

You can immediately see that the maximum value rose above 1.0 and mass changed. Those are useful clues.

Logging vs Print Statements

For serious debugging, use Python’s logging module instead of raw print().

import logging

logging.basicConfig(
    level=logging.DEBUG,
    format="%(levelname)s: %(message)s"
)

logger = logging.getLogger(__name__)

logger.debug("Running diffusion step at t=%f, dt=%f", t, dt)
logger.warning("Energy increased by %.5f", energy_change)

Logging has several advantages:

  • Levels such as DEBUG, WARNING, and ERROR.
  • Optional timestamps and file output.
  • Easy activation or deactivation without deleting debug code.
  • Cleaner diagnostics during long simulation runs.

A practical rule is to use print() for quick exploratory scripts and logging for reusable or production-style simulation code.

Level 2: Assertions and Sanity Checks

Assertions are like print statements that refuse to stay silent. They raise an error when an assumption is violated.

The Pattern: Assumptions First

The most useful pattern is to check inputs, run the solver, and check outputs.

import numpy as np

def solve_diffusion(phi, dt, diffusion_coeff):
    # Preconditions
    assert np.all(phi >= 0), f"Negative phase: {np.min(phi)}"
    assert np.all(phi <= 1), f"Phase > 1: {np.max(phi)}"
    assert np.isfinite(phi).all(), "Non-finite input values detected"
    
    # Solver step
    phi_new = apply_diffusion(phi, dt, diffusion_coeff)
    
    # Postconditions
    assert np.isfinite(phi_new).all(), "Solver produced NaN or Inf"
    assert np.all(phi_new >= 0), f"Negative output: {np.min(phi_new)}"
    
    return phi_new

This turns silent numerical failure into a clear error with location and context.

The Assertion Philosophy

Assertions are not only about catching bugs. They pin bugs down.

An assertion that fails near the start of a solver tells you where the assumption broke. A NaN that appears hundreds of lines later tells you much less.

Useful places for assertions include:

  • Before and after solver calls.
  • At module boundaries where data is exchanged.
  • After mesh operations such as refinement, remeshing, or coarsening.
  • After applying boundary conditions.

A Warning About Assertions

Do not overuse assertions in expensive inner loops unless you need them. They can add runtime overhead and make debugging noisy.

Use assertions for invariants and preconditions: conditions that should never be violated if the code is correct.

Level 3: The Python Debugger

When print statements are not enough, use an interactive debugger. A debugger lets you pause the program, inspect state, and step through logic.

The Classic Tool: pdb

pdb is Python’s built-in debugger. It is available in every Python installation.

In IPython, you can use post-mortem debugging after an error:

In [1]: %run simulation.py

In [2]: %debug

This drops you into the debugger at the point of failure.

For command-line scripts, launch the script with the debugger:

python -m pdb simulation.py

Useful commands include:

  • c: continue execution.
  • n: move to the next line.
  • s: step into a function call.
  • b: set a breakpoint.
  • p variable_name: print a variable.
  • q: quit the debugger.

You can also insert a breakpoint directly in code:

def solve_phase_field(phi, dt):
    import pdb
    pdb.set_trace()
    
    phi_new = update_phase_field(phi, dt)
    return phi_new

ipdb and IPython Debugging

ipdb is an enhanced version of pdb. It provides tab completion, better formatting, and a more comfortable interactive experience.

pip install ipdb

import ipdb
ipdb.set_trace()

For IPython users, %debug is often the fastest option after an exception.

IDE Debuggers

PyCharm, VS Code, and Spyder include graphical debuggers. They let you set breakpoints, inspect variables, step through code, and watch expressions without typing commands.

IDE debuggers are helpful for local development. For remote clusters or headless servers, pdb, ipdb, and logging are usually more practical.

When to Use Which Tool

Scenario Recommended Tool
Quick script or early local development print() plus assertions
Error traceback in IPython %debug
Remote cluster without IDE access python -m pdb script.py
Interactive debugging with better context ipdb or IDE debugger
Production-style diagnostics logging with strategic breakpoints

Level 4: Profiling to Find Hidden Bugs

Profiling is usually seen as a performance tool, but it also helps debugging. Slow code is often wrong code. Profilers can reveal accidental Python loops, repeated allocations, redundant calculations, or unexpected function calls.

cProfile: Find the Hot Functions

import cProfile
import pstats

profiler = cProfile.Profile()

profiler.enable()
solve_phase_field(phi, dt, n_iterations)
profiler.disable()

stats = pstats.Stats(profiler)
stats.sort_stats("cumulative")
stats.print_stats(20)

This shows which functions consume the most total time. If a helper function unexpectedly dominates runtime, it may contain a hidden bug or inefficient array logic.

line_profiler: Find the Hot Lines

When cProfile identifies a slow function, line_profiler can show which lines inside the function consume time.

pip install line_profiler
kernprof -l -v simulation.py

This is useful when a single line inside a loop creates arrays repeatedly or performs a slow operation by mistake.

Memory Profiling: Find Hidden Allocations

Scientific code often hides memory bugs such as accidental copies or excessive array creation inside loops. memory_profiler can help track allocations.

from memory_profiler import profile

@profile
def solve_phase_field(phi, dt):
    phi_new = apply_diffusion(phi, dt)
    return phi_new

Line-by-line memory profiles can reveal where large arrays are created and discarded.

Profiling as Debugging

Profiling helps when a solver is both wrong and slow. For example, a wrong slice may compute energy over the wrong region. A profile can show that the energy calculation is unexpectedly expensive, which points you back to the array operation.

Assertions confirm the bug. Profiling helps you find where to look.

Level 5: Testing and Regression Prevention

The final debugging level is prevention. After fixing a bug, write a test that reproduces the failure and proves the fix works.

The Scientific Testing Pattern

import pytest
import numpy as np

@pytest.mark.parametrize("mesh_size", [10, 20, 50])
def test_energy_decreases(mesh_size):
    """Energy should decrease during simple diffusion."""
    phi = init_phi(mesh_size)
    energy_history = []
    
    for step in range(100):
        energy_history.append(compute_energy(phi))
        phi = solve_step(phi, dt=0.01)
    
    energy_changes = np.diff(energy_history)
    assert np.all(energy_changes <= 1e-10)

This test checks physical behavior, not only syntax. A correct simulation should not merely run. It should conserve expected quantities, respect bounds, and converge with refinement.

The Testing Pyramid for Scientific Code

  • Unit tests are fast and check individual functions such as compute_gradient() or apply_boundary_conditions().
  • Integration tests run solver-level checks against analytical or manufactured solutions.
  • Regression tests compare full simulations against trusted baseline outputs.

The key principle is simple: test the physics, not only the code path.

Common Scientific Debugging Pitfalls

NaN and Infinity

NaN and infinity values are common in scientific code. They often appear after division by near-zero values, log(0), sqrt of negative values, or overflow in exponentials.

phi = np.zeros((100, 100))
phi[50, 50] = 1.0

result = phi / (phi + 1e-300)

Check for finite values after numerical operations:

assert np.isfinite(result).all(), (
    f"NaN or Inf detected: min={np.nanmin(result)}, max={np.nanmax(result)}"
)

Array Shape Mismatches

Array shape errors are common when moving between flattened solver arrays and 2D or 3D fields.

print(f"phi shape: {phi.shape}")
print(f"expected shape: {expected_shape}")

Common causes include accidental flattening, indexing that removes a dimension, and library variables that store data in a different shape than expected.

Floating-Point Precision

Do not compare floating-point arrays with exact equality unless exact equality is truly expected. Use tolerance-based comparisons.

np.testing.assert_allclose(phi, phi_new, rtol=1e-6, atol=1e-12)

This is especially important for iterative solvers, mixed precision, and long time integrations.

Boundary Condition Bugs

Boundary condition bugs are common because the code often runs even when the wrong cells or faces are constrained.

# Wrong: overwrites all cells
phi[:] = boundary_value

# Better: apply only to selected boundary locations
phi[boundary_mask] = boundary_value

Always verify that boundary values are applied only where intended and that interior values are not overwritten.

Convergence and Stability

Simulations can diverge silently when time steps are too large or solver tolerances are too loose. Monitor conserved or bounded quantities over time.

if energy > prev_energy * 1.1:
    print(f"Energy growing. Reducing dt from {dt} to {dt / 2}")
    dt /= 2

Track energy, mass, momentum, bounds, residuals, and time-step behavior. Plot them during debugging.

The Debugging Escalation Ladder

Use the lowest debugging level that gives enough information.

Level 1: print()      — quick inspection
Level 2: assert()     — sanity checks and invariants
Level 3: pdb/ipdb     — interactive state inspection
Level 4: profiling    — performance bugs and hidden inefficiency
Level 5: pytest       — regression prevention

Escalate when needed:

  • If print() gives enough context, stay at Level 1.
  • If you need to verify invariants, use assertions.
  • If you need to inspect state at a specific line, use pdb or ipdb.
  • If the code is wrong and slow, profile it.
  • If you want to prevent the bug from returning, write a test.

A Practical Debugging Example: Fixing a FiPy Simulation

Suppose a FiPy diffusion simulation produces unbounded concentration values instead of smooth diffusion. Start by writing a minimal, clear version of the problem.

import numpy as np
from fipy import Grid2D, CellVariable, TransientTerm, DiffusionTerm

mesh = Grid2D(nx=50, ny=50, dx=1.0, dy=1.0)

phi = CellVariable(name="concentration", mesh=mesh, value=0.0)
x, y = mesh.cellCenters

# Initial hot spot in the center
phi.setValue(
    1.0,
    where=((x - 25.0)**2 + (y - 25.0)**2) < 25.0
)

eq = TransientTerm(var=phi) == DiffusionTerm(coeff=1.0, var=phi)

dt = 0.01

print(f"Initial: min={phi.value.min()}, max={phi.value.max()}")

for step in range(100):
    eq.solve(var=phi, dt=dt)
    
    assert np.isfinite(phi.value).all(), "NaN or Inf detected"
    assert phi.value.min() >= -1e-10, f"Negative value at step {step}"
    
print(f"Final: min={phi.value.min()}, max={phi.value.max()}")

Then debug systematically.

Step 1: Make It Fail Reliably

Log minimum, maximum, mass, and energy before and after the run.

print(f"min={phi.value.min()}, max={phi.value.max()}")
print(f"mass={np.sum(phi.value * mesh.cellVolumes)}")
print(f"energy={np.sum(phi.value**2 * mesh.cellVolumes)}")

Step 2: Divide and Conquer

Test the diffusion term alone. Remove reaction terms, advection terms, nonlinear sources, and complex boundary conditions until the failure disappears or becomes isolated.

Step 3: Change One Thing

Try a smaller time step, then compare results.

dt = 0.001

If the bug persists with a smaller time step, the problem may be spatial, structural, or related to boundary conditions rather than time stability.

Step 4: Use the Debugger

Drop into a debugger near the failing step and inspect values, mesh shape, masks, and boundary conditions.

import pdb
pdb.set_trace()

Step 5: Add a Test

After fixing the bug, add a regression test.

def test_diffusion_mass_stays_finite():
    phi_initial, phi_final, volumes = run_small_diffusion_case()
    
    mass_initial = np.sum(phi_initial * volumes)
    mass_final = np.sum(phi_final * volumes)
    
    np.testing.assert_allclose(mass_initial, mass_final, rtol=1e-6)
    assert np.isfinite(phi_final).all()

Summary and Next Steps

Debugging scientific Python code is a progressive skill. Start with print statements, move to assertions, use debuggers when state inspection is needed, profile when the code is slow or suspicious, and write tests to prevent regressions.

The most important principle is to isolate the bug before fixing it. A bug that you can reproduce in a small test case is much easier to solve. A bug that only appears inside a large simulation remains a mystery.

For scientific code, remember these rules:

  • Check for NaN and infinity after numerical operations.
  • Print array shapes when indexing becomes complex.
  • Monitor conserved quantities such as mass, energy, and momentum.
  • Test boundary conditions separately from interior logic.
  • Use assertions as the first line of defense against silent numerical failure.

Next Steps

  • Start using the debugging escalation ladder: begin with print() and escalate only when needed.
  • Add assertions before and after solver functions.
  • Set up logging for simulation diagnostics.
  • Write at least one test case that reproduces a common bug.

Debugging is not a sign of weakness. It is part of writing code that matters. The skills you build while debugging scientific Python code will help across every kind of numerical and research software project.

Related Guides

Need Help Debugging Your Simulation Code?

Debugging complex scientific simulations can be time-consuming, especially when you are unsure whether the bug is in the code, mesh, boundary conditions, or solver settings.

Our team can help with systematic debugging workflows, numerical test cases, sanity checks for conservation laws, and CI pipelines with regression tests.