Why Your Simulations Keep Exploding
The first time I tried to model a simple pendulum with drag, my solution diverged so fast it took three seconds to crash. A full hour went into debugging before I realized the timestep was way too big for the friction term. That kind of thing is why I keep these tricks around. They are not glamorous. They just keep the numbers from becoming nonsense. I do not write a single simulation without checking a few things before I let it run. The process is boring because the payoff is that I stop wasting time staring at NaNs. I start with the units, move to the timestep, test the boundary conditions, and then add damping only where it makes sense. Here is the order I use. First, I verify the model by nondimensionalizing the governing equations. That is not just academic padding. When I write the drag coefficient and the Reynolds number into the same block, I can see exactly which term is responsible for the blow-up. Second, I set the timestep from the fastest physical timescale. If the fastest mode is something like the acoustic transit time across a single mesh cell, the timestep needs to be small enough to resolve it. Third, I run a zero-length test with no source terms and no boundaries. The system should be static or moving at constant speed. If it is not, the initial conditions or the discretization are wrong.
I learned this the hard way with a heat conduction problem that used an explicit scheme. The stability limit is alpha times dt over dx squared. I thought my mesh was fine. It was not. I got a nice-looking temperature profile until the eighth step, when it went wild. I reduced dt by a factor of ten and the profile looked normal. The workaround is always the same. Use the explicit stability limit as a ceiling, then drop it by another factor to give yourself margin. I also make sure my boundary conditions are consistent with the physics. A fixed temperature at one end and a zero flux at the other is fine. A fixed temperature at both ends with a source term in between is also fine. A fixed temperature at one end and a flow rate at the other for an incompressible fluid is a different story. You need to think about what information the solver needs and what it is allowed to compute. Otherwise you get an underdetermined or overdetermined system. Both fail in different ways. Energy conservation is the single most useful sanity check. If I integrate the power input over time and compare it to the change in kinetic plus potential energy, the numbers should match within a small tolerance. In practice they do not match perfectly because of numerical dissipation, but they should be in the same ballpark. I once modeled a spring-mass system and the energy grew over time. That told me the integrator was adding energy. I switched from a fourth-order Runge-Kutta to a symplectic method and the drift stopped. That is a common fix for long-running orbital and vibration problems.
I use a handful of standard tools to make this reproducible. Python with SciPy does most of the heavy lifting. I write the equations in a vector form and let the integrator handle the time stepping. If I need more speed, I move the core loop to C or use Numba. For larger meshes, I use deal.II or FEniCS. For particle systems, I use LAMMPS. None of this is novel. It is just what works for me.
Get the Full Details
A Real Problem I Faced and the Workaround
Last year I worked on a two-phase flow problem with a dense particle phase. The drag term was stiff because the particle response time was much smaller than the fluid timescale. I tried a standard explicit coupling. It failed after about five iterations. The solver kept oscillating and the density field went negative. I switched to an implicit coupling with a subcycling approach. I solved the fluid equations on the main timestep and updated the particle velocity inside each fluid step with a smaller internal step. The code got longer. The results stopped blowing up. That is the pattern. Stiff coupling usually needs an implicit treatment. I also ran into a problem with contact forces in a rigid-body simulation. The impulse-based method produced energy spikes when two objects collided at a shallow angle. I added a small restitution coefficient and used a penalty method with a softened contact stiffness. The spike dropped to a manageable level. This is not a perfect fix. It changes the physics slightly. But it is often the only practical fix when you cannot afford a full DEM solver.
When These Shortcuts Fail
I need to be honest about the limits. These tricks work well for low-Mach-number flows, linear elasticity, and moderate Reynolds numbers. They break down when you have strong shocks, turbulent cascades at high Reynolds numbers, or highly nonlinear material behavior. A good example is supersonic flow past a blunt body. The explicit methods become unstable because the Mach number is high. You need an implicit scheme or a Godunov-type solver with Riemann solvers. Another example is granular flow with large strain. The standard continuum models smear the failure surface. You need a micropolar or peridynamic model. Those are much more expensive and harder to tune. There is also a practical limit to how much you can simplify the mesh. If the geometry has small features, like a thin boundary layer, you need many cells there. A coarse mesh will give you wrong results even if the solver is stable. I usually refine the mesh until the solution stops changing. That is the only reliable check. Analytical solutions are rare in real problems. Numerical experiments are the norm. If you are just starting out, I recommend you begin with 1D problems. The wave equation in one dimension is simple enough to solve by hand and easy to implement. Then move to 2D heat conduction. Then try a 2D incompressible flow with a known solution. Once you can reproduce the known solution, you can trust the unknown ones. That is the only path I have found that does not waste time.
Resources and Where to Find Code
I do not maintain a single repository. I keep examples scattered across my personal GitHub and a few public repos. If you want to follow along, look for simple solvers that implement the ideas above. The code is not polished. It is meant to be read and modified. I include comments that explain why I chose a particular discretization. That is more useful than a clean API. For learning, I recommend you read the basics of numerical analysis. Not all of it. Just the chapters on stability, consistency, and convergence. Then read a book on finite elements or finite volumes. A good starting point is the work by Zienkiewicz or the classic papers on stabilized methods. After that, you can pick a library and start breaking things. That is how you learn what fails and why. I also keep a notebook of edge cases. One entry describes a situation where the pressure Poisson equation became ill-conditioned because the domain had a floating body. The fix was to add a small pressure constraint at a single point. Another entry describes a case where the mass matrix was singular because a node was unconstrained. The fix was to add a dummy spring. These are the kinds of problems that do not show up in tutorials.
If you want to experiment with these concepts, you can download the example code I mentioned. It is not a product. It is a collection of scripts that demonstrate the tricks. You will need to adapt them to your own problems. That is the point. You should not copy-paste. You should understand. The scripts use standard libraries. No special compilers. No paid licenses. Just what you can run on a laptop. The main takeaway is simple. Stability is not optional. Conservation is a useful check. Mesh refinement is necessary but not sufficient. And the easiest way to learn is to break things deliberately and see what happens. I have spent years doing exactly that. The shortcuts I described are the ones that survived the breaking. I also keep a list of common pitfalls. One is assuming that a smaller timestep always improves accuracy. It does not. In an explicit method, it can improve stability. In an implicit method, it can increase roundoff error. The optimal timestep depends on the problem. Another pitfall is ignoring the boundary condition implementation. A wrong boundary condition can dominate the solution even if the interior scheme is perfect. I check the boundary conditions by comparing the numerical result to an analytical solution at the boundary. If they do not match, I revise the implementation.
Finally, I want to mention that these tricks are not a substitute for understanding the underlying physics. They are tools to make the simulation work. The physics still drives the result. If the physics is wrong, the tricks will not save you. You need to know what you are modeling. That is the part that takes the most time. The rest is just coding and testing.
