The Coastline Broke My Electromagnetic Solver, So I Rebuilt It
The last time I wrote about this, I thought I’d finally figured out why coastlines and the seafloor were breaking my electromagnetic (em) solver. The fields are represented with spherical harmonics, which are really good at smooth things. But, the Earth is... not smooth. There’s water here, rock there, and at these boundaries, conductivity suddenly changes.
So, instead of throwing more and more numerical points at those sharp boundaries, I decided to actually find the boundary inside each computational cell, split the cell, and integrate the water and rock separately. And, I really thought that might finally solve the convergence problem.. but, no it didn’t :(
What am I actually solving?
The physics itself is easy: A changing magnetic field creates an electric field, and that electric field drives current through conductive material. Some of that electrical energy gets dissipated as heat.
∇ × E = iωμH
J = σE
q = ½σ|E|²
The hard part is convincing my computer to solve it across an entire spherical ocean + solid Earth without creating numerical garbage every time water turns into rock. This is what I do for fun now. :)
G: handle the boundary better
Phase 3G was a series of attempts to integrate across sharp conductivity boundaries more accurately. I was trying to see if a computational cell contains both seawater and rock, can I still calculate the material interaction accurately enough without rebuilding the entire solver.
I tried heavier dealiasing, and a much finer quadrature. None of it worked consistently enough, and so at that point I had a choice. I could keep dialing up the resolution until the result looked stable enough...or admit the numerical treatment still wasn’t good enough.
So, fine, it’s not good enough. Eventually I want to compare this heating pattern to the actual abyssal ocean, so it has to hold up. So, I moved on to Phase 3H..
H: actually represent the boundary
This was the idea from the end of my last post. I made the integration conform to the actual interface, or in other words, I found the boundary and integrated each side. To do this, I built a cut cell treatment: if a cell crosses the water-rock interface, I split it. Water gets integrated as water, and rock gets integrated as rock. The conductivity of each remains separate, and not blurred.
And, this part worked.. the calculation finally became accurate enough to move on. Of course then it became computationally ridiculous. 16,000 years to complete, so we’d get an answer in like almost 3 geophysical events :)
So, now we have an engineering problem: how do I keep the exact interface treatment without constructing enormous matrices and asking too much of my machine?
I went through a few versions and eventually landed on a row wise construction that fits into the actual production workflow. So, in 3H, we definitely fixed the material boundary problem, but now we have a convergence problem.
I: please converge
Even after fixing the interface calculation, the full electromagnetic solution still wasn’t converging the way I wanted. The total energy looked pretty good, but the detailed electric field did not. Total bummer because I care about where the energy goes, especially at depth.
I decided I needed an independent benchmark. Basically, I needed another solver to check the first solver. So, I built an H(curl) finite element benchmark. Finite elements chop the volume into a huge number of little 3D pieces, and in our case: little baby tetrahedra. The important part is that I can make the water rock boundary lie directly along the faces of those tetrahedra. So the material interface is not something the solver has to infer from samples, it’s literally built into the mesh.
Now I can actually determine if the spherical solver is wrong, or if this is just a genuinely difficult convergence problem. I’m using a fake Earth for this test, nothing fancy yet. Can we just show this works?
Earth radius = 6,371,200 m
Modeled depth = 12 km
Seawater conductivity = 3.2 S/m
Solid conductivity = 0.01 S/m
So, the finite element mesh started behaving (excellent), but the electromagnetic boundary condition on the inner and outer surfaces of the shell also uses spherical harmonics.
So I fixed one boundary... and immediately had to qualify another one.
There’s another harmonic cutoff here, called L. Higher L means the electromagnetic boundary can represent finer angular structure. At first I could only support a fairly low cutoff (not enough), so I increased it, then increased it again.
larger L = more angular detail available at the boundary
Eventually I had to redesign the whole boundary trace system to get where I needed to go.
And, yeah, I though about quitting…. until the new boundary representation succeeded. It can now represent every spherical harmonic mode I need through degree 108! And, this time the numerical tests look really good.
If I give the system a known harmonic mode, project it onto the finite element boundary, and ask it to recover that mode again...basically floating point noise (excellent).
Originally, I had also thought connecting all of the finite element boundary degrees of freedom to all those spherical harmonic modes might require one enormous dense matrix. Absolutely did not want that. Fortunately, the spherical symmetry lets me separate the problem into independent m blocks. So instead of one gigantic projection, I solve lots of small pieces.
Surprisingly, it’s fast. Both the inner and outer boundary operators passed linearity, repeatability, and adjoint tests essentially to machine precision. So the boundary calculation is not the bottleneck anymore (huge win).
Once the L108 boundary passed, I could finally build the full finite element benchmark mesh that supports it. And, this is where I looked at the numbers and thought...ohhhh 😂
More than a million electromagnetic unknowns, great. But, the mesh passed all of its checks.
where I am now
So after all of that, I finally ran the full electromagnetic problem from L32 all the way through L108. And... it still didn’t converge 😂 But, I still learned something.
The total Joule heating is basically stable now. From L32 to L108, the total power only changed by about 0.08%, and every step passed the power-convergence test. This means the total watts are not really the problem anymore. It’s the electric field that’s the problem.
At L108, the field is still changing by about 0.00215 compared with the previous run, and my target is 0.0005, so I’m still a little more than 4x above where I want to be.
And, this is where I decided to stop just increasing the resolution, because I could keep pushing L higher and higher. But, I don’t want to just throw more resolution at the solver until eventually the number gets small enough.
I want to know why it’s still changing. So, I went back into the saved solutions and looked at where the remaining error is actually showing up, and I think you’ll find this interesting.
At the boundary, about 86% of the remaining correction is still sitting in the newest spherical harmonic degrees I’m adding. This means increasing L is definitely still changing something in the boundary representation.
But deeper inside the Earth, the pattern looks different. At the highest resolutions, most of the remaining electric field change is more radial and increasingly concentrated deeper in the 12 km shell. In the final L106 → L108 step, about 86% of the remaining field change energy is between 4 and 12 km depth.
So now I have two possible paths: 1) I still need more spherical-harmonic detail at the boundary, and 2) the finite element mesh inside the Earth has reached its own resolution limit, and I’m starting to see that instead.
And with only one mesh, I can’t tell which one it is. So, that’s the next test. I’m building the exact same L108 problem again, but this time with finer resolution through the depth of the Earth. Then, I compare the two.
If the electric field barely changes, that tells me the mesh was probably already fine enough and I have a much stronger reason to extend the harmonic boundary beyond L108. If the field changes a lot, especially in that same deeper radial pattern, then I’ve probably found the actual limit of the finite element mesh.
And in that case, increasing L would’ve just meant throwing more spherical harmonics at the wrong problem. Either way, I get an answer.
Which, after rebuilding basically every numerical part of this solver at least once, is what progress looks like :)
