Stop Building Spheres Out of Bricks
How a master’s student cut the error of an acoustic simulation by more than 300 times at the finest grid, just by changing the elements next to a curved wall
There is a number in my student’s results table that I keep coming back to: 274.6 percent.
That is the error of a simulated echo. An error above 100 percent means the simulation is off by more than the signal it is supposed to compute. It is not “a bit inaccurate”. It is a different echo.
Then he changed how the solver treats the cells around a curved wall, ran the same case on the same grid, and the error dropped to 5.0 percent. On the finest grid we tested, it went from 68.9 percent to 0.21 percent.
In this blog, I want to walk you through what happened. It is a story about a master’s thesis, about spheres, and about a lesson every simulation person learns eventually: you cannot refine your way out of a wrong shape.
What the thesis is trying to do
I supervise a master’s student, a colleague from work, whose thesis deals with ultrasonic parking sensors. These are the little discs in a car bumper that beep when you approach a wall. They work at roughly 40 to 50 kHz, and companies want to train machine-learning models on their raw echo signals. Those models need an enormous amount of labelled data, far more than anyone can economically record on real roads. So the data has to come from simulation.
There are two ways to simulate sound, and neither is enough on its own.
Wave solvers solve the acoustic equations directly. Diffraction and interference simply fall out of the physics. The price: the thesis task puts a single 3 × 2 × 1 m scene at 40 kHz at about 76 billion unknowns and roughly 25 CPU-years. You do not generate thousands of scenes like that.
Ray tracing is fast. It sends rays through the scene, like light. But it is an approximation, and it is least reliable exactly where a parking scene gets interesting: sharp edges, narrow gaps, and surfaces curved at about the size of the wavelength.
So the thesis builds a hybrid. Solve the wave problem once per canonical object (a sphere, say) with an accurate wave solver, and hand the result to the ray tracer. Part of the work is then to find out where this hybrid can be trusted.
That plan only works if the wave solver is right in the first place. The first research question in the thesis is therefore almost boring in its honesty: how fine do the elements and the polynomials have to be before the echo stops changing?
And that is where the elements come in.
Discontinuous Galerkin in two minutes
The wave solver uses the Discontinuous Galerkin method, or DG. I worked with it during my PhD, so it is close to my heart.
The idea: chop the domain into small elements. Inside each element, describe the pressure and the velocity with a polynomial. Unlike the classic Finite Element Method, the polynomials of neighbouring elements are not forced to match at the edges. Instead, the elements talk to each other through a flux across each shared face, an idea borrowed from the Finite Volume method.
Two knobs control accuracy:
- h, the size of the elements
- N, the polynomial degree inside each element
Smaller h or larger N should give a better answer. For smooth problems, degree-N elements are expected to converge at roughly order N+1: halve h, and the error drops by about 2^(N+1). For N = 3, that is a factor of 16 per halving.
Should is doing a lot of work in that sentence.
The problem: spheres made of bricks
Our elements are cubes. A sphere is not made of cubes.
The original solver handled this the simplest way possible: any cube whose centre lies inside the sphere is treated as solid, every other cube as air, and the sphere becomes a staircase of cubes. Think of a sphere built from Lego bricks. We call it the “staircase” approach.
At our coarsest grid, the sphere has a radius of 20 mm and the elements are 10 mm wide. By my arithmetic from the recorded degrees of freedom, the staircase sphere is a blob of 32 cubes. To a sound wave, that is a very different object. Instead of one smooth reflection, you get reflections off flat faces and corners, a different wall position, and different diffraction.
Here is the sneaky part. The polynomial inside each cube can be as accurate as you like. It does not matter. The wall is in the wrong place. That is a geometry error, and a staircase’s geometry error shrinks only about linearly with h. It puts a ceiling on how good the answer can get, however high you push N. The data show this nicely. At a 2⁴³ grid, going from N = 2 to N = 3 made the staircase result no better (96.8% to 186.5% error). With cut cells, the same step took it from 17.7% to 0.74%.
You want to read more blogs? Become a member and support me and thousands of other writers directly. If you use my link I get paid a monthly commission, however, at no extra cost for you.
The fix: let the sphere cut through the cells
The alternative is called a cut-cell method. Instead of rounding the sphere to the grid, you let the sphere slice through the cubes and keep the part of each cube that is still air. The wall is now exactly where the sphere is, with the exact surface normal at every point.
In my student’s solver, that means:
- The integrals over a cut cell run only over its air portion, using quadrature that follows the true sphere.
- The wall condition uses the analytical sphere normal, not the normal of a brick face.
- Nothing gets filled in or thrown away: every sliver of air counts.
This is easy to say and nasty to do. Cutting a cube with a sphere can leave an absurdly thin sliver of air, and thin slivers wreck explicit time stepping. The cell has almost no volume, so its equations become extremely stiff and the stable time step collapses. This is the well-known “small cell problem”. At the finest grid, the smallest air fragment is 0.188 percent of a full cell.
Get Beyond Boundaries’s stories in your inbox
Join Medium for free to get updates from this writer.
The standard escape is aggregation. Every cut cell shares one polynomial with a neighbouring full air cell (its “host”). Picture a tiny sliver of cake that you glue onto a whole slice, so that the two are treated as one piece. The sliver no longer carries its own fragile unknowns, and the stiffness problem goes away without throwing any air away. The solver was also built so that the discrete energy cannot grow at the wall, which is exactly what you want from a rigid boundary.
What the results show
Time to see whether all of this pays off.
The test: a 16 cm cube of air, a rigid sphere of 20 mm radius, a source and a receiver close together, and a short ultrasonic pulse. The solver runs twice, once with the sphere and once without. Subtracting the two leaves only the echo from the sphere. That echo is then compared against a saved semi-analytical reference. Nothing is rescaled, shifted, or fitted. The metric is the relative RMS difference of the echo envelope, N = 3 throughout:
| Grid (element size) | Staircase | Cut cells | Cut cell is better by |
| 16³ (10 mm) | 274.6% | 5.01% | about 55 times |
| 24³ (6.7 mm) | 186.5% | 0.74% | about 250 times |
| 32³ (5 mm) | 68.9% | 0.21% | about 330 times |
Three things stand out.
The staircase is not just inaccurate, it is erratic. On the coarsest grid its echo peaks at 0.069 Pa, about 2.5 times the reference peak of 0.0280 Pa. On the finest, its echo splits into two almost equal bumps, about 0.2376 and 0.3036 ms, with an amplitude ratio of 0.93. Its error keeps falling, but the solutions barely agree with each other: the difference between the 24³ and 32³ staircase envelopes (146%) is larger than between 16³ and 24³ (56%). When successive refinements do not settle down, the discretisation is not converging, no matter what the error column says.
The cut-cell solutions settle down. The same two numbers go from 4.37% to 0.76%. The peak amplitude sits at 0.0262, 0.0278, and 0.0281 Pa against a reference of 0.0280 Pa. The peak time matches the reference to within one sample of the reference time axis at every grid. (I would not call that sub-sample timing accuracy, and neither does my student.)
The local slopes look like what theory suggests. My back-of-the-envelope fit of the cut-cell errors against h gives local slopes of roughly 4.7 (16³ to 24³) and 4.4 (24³ to 32³). For degree-3 elements, N+1 = 4 is the ballpark. But three grids, with a reference that has its own model error, do not make a convergence order. The report makes no such claim, and I won’t either.
What did it cost?
The honest answer: some. The cut-cell runs on the same cell count are slower (sphere runtimes: 146 s vs 70 s at 16³, 474 s vs 347 s at 24³, 1,231 s vs 1,052 s at 32³). Both methods have almost the same number of unknowns per grid, 258,048 vs 260,096 at 16³. These timings were taken while other jobs were running, so they are not controlled benchmarks.
But here is the comparison I like best. The coarse 16³ cut-cell run took about 146 seconds and had 5.0% error, with roughly 258,000 unknowns per field. The finest staircase run, 32³, took about 1,052 seconds, had about 2.08 million unknowns per field, and showed 68.9% error. Eight times more unknowns, seven times the runtime, and a much worse echo.
Can we trust it?
A fair question, and my student spent a lot of time on it. In the final report:
- In a closed domain without forcing, the energy never exceeded its starting value in the tests. The global pressure conservation error was below 1.3e-17.
- Halving the time step changed the echo by 0.0086%.
- Raising the quadrature order for the cut geometry changed it by about 5e-6%.
- A flat-wall check against the method of images gave a 0.0002% waveform difference in the tested configuration.
And the limits, stated plainly:
- This echo is dominated by about 11.6 kHz, not 40 kHz. The Gaussian source filters the carrier so strongly that the 40 kHz content is minuscule. The agreement shows the geometry treatment works. It does not yet prove 40 kHz resolution, which is what the thesis ultimately needs.
- The reference is itself a model (an approximate source and an unbounded smooth sphere), so part of the remaining difference is model difference.
- Aggregation changes the polynomial space near the wall too, so this is not a pure geometry-only comparison.
- A shifted sphere needs higher quadrature, and an optional finer grid and N = 4 pulse runs were not done.
- It is a CPU-only experimental solver, not a validated tool for arbitrary shapes.
What I take from this
The headline is not “cut cells are magic”. It is that the shape of your discretised wall is part of your numerical method, and sometimes the dominant part. We were refining h and N on a sphere that was never really a sphere. Every refinement was spent polishing the wrong object.
If you work with simulations, here is a practical check you can run before burning more compute:
- Refine the mesh and see whether the answers get closer to each other, not only to a reference.
- Raise the polynomial degree and see whether anything changes.
- If neither helps much, suspect the geometry first.
My student’s next steps are the unglamorous ones: 40 kHz content, more shapes, and eventually the hybrid with the ray tracer and the lab measurements. But the foundation now stands on a sphere that is actually round, and I find that quietly satisfying.
Remember the 274.6 percent? It was never a bug in the physics. It was a pile of bricks pretending to be a sphere.
Thank you very much for reading. If you liked what you read, I would appreciate a clap 👏, and don’t forget to follow me on Medium and Twitter for more blog posts.
Would you like to receive my free weekly newsletter „The Creative Engineer”? Make sure to subscribe here.
You have feedback, complaints, want to collaborate, or anything you want to ask me? Please don’t hesitate to write me a DM on Twitter or an e-mail at v.ganiu@gmail.com. Or just leave a comment on this blog.
Would you like to support me buying an USB-C hub for my desk setup? You can buy me a coffee or donate to my PayPal. Thank you very much for any kind of support!
