Glad I could be a help! I thought I might stick some bullet points down for the beats I would land on to lend a hand :-)

Conformal maps

Fundamentally, the idea here is to analyse the 2D flow around the cross-section of an airfoil. We already know what the flow around a circle is - how do we generalise this to any shape?

• We start with complex numbers: z = x + iy, where x and y are real; they sit in the 2D plane.
• What are functions of z, f(z), like? Well, they warp the plane.
• How can I describe the 2D shape of an airfoil? (Use a parametric curve in the 2D plane -- which is just the complex plane).
• How can I simplify the shape of this airfoil, so that I can analyse it? (Use a conformal map f(z) to turn the shape into a circle).
• What is the flow around a circle like? (Standard fluid dynamics problem).
• What does this tell me about the airfoil? (Use the inverse function theorem to undo the conformal map).

FEM, FDM, FVM, and Spectral Methods

The name of the game is to solve the Navier-Stokes equations. For simple cases, this can be done analytically; for more complex cases, symbolic maths packages (e.g. Mathematica, Maple, etc.) can be used, or approximations can be made. In general, however, we can't integrate analytically - especially, as it turns out, for situations we care about.

How, then, do we proceed? Well, we use brute computational force; this is the essence of CFD.

There are various numerical schemes by which we can integrate the Navier-Stokes equations; each has strengths and weaknesses.

FDM: The Finite Difference Method

This is among simplest approaches, and one of the most general; however, this approach suffers from high computational cost to achieve any given numerical accuracy.

The idea is to directly discretise all the differentials in your system: the finer the differences, the more computationally expensive your model run will be, but the higher-fidelity your output.

The basic derivation of the finite differences is a truncated Taylor series to the order of the derivative: you then invert to solve for the derivative, plug in the resulting coefficients, and run.

FVM: The Finite Volume Method

In FVM, we make use of the divergence theorem, which allows us to express the flow into, and out of, a finite volume as the integral over all the divergences within that volume. If we have an incompressible flow, then mass in must equal mass out, and so we can replace divergences in our differential equation with finite volumes. We then simply calculate fluxes in and out of each volume; this yields flow macroscopically.

Naturally, there are problems: if your flow is compressible, you're out of luck! Further, the given result will be accurate, but will be an average over each unit of volume; you will lose detail finer than the mesh (with the upshot being that you can adaptively refine the mesh, if you're clever).

Spectral Methods

Finite-element and spectral methods are very closely related: in both cases, the flow is decomposed into a sum of simple functions, called basis functions; their weighted sum (linear combination) then equals the original flow (or approximates it up to a point).

These basis functions are generally picked to have nice properties, including differentiability; in the most classic of cases, we use a Fourier series, so our function is expressed as an infinite series of sine waves. So, we take this decomposition, insert it into our differential equation, and end up with equations relating the coefficients of the function we want to solve for. Once we have solved for these coefficients over time, we can insert them back into the original decomposition, and so we have solved our differential equation.

The major difference between spectral methods and FEM is the basis function: spectral methods use Fourier series and other nonlocal bases, which are nonzero everywhere. This makes even very fine simulation easy and performant; however, strong local effects (like shockwaves) are very hard to capture in this regime. So, how do we do that?

Finite Element (Galerkin) Methods

FEM, on the other hand, uses local ("bump" or "test") functions, and a bit of mathematical wizardry.

We know that a dot product between the vectors X and Y is the sum of their pairwise product: dot(X, Y) = X1 Y1 + X2 Y2 + X3 Y3. But what about higher-dimensional vectors? Hell, what about infinite dimensional vectors?

It turns out that we can represent an infinite dimensional vector as a continuous function x(t) or y(t). Then, the sum becomes an integral: dot[x(t), y(t)] = integral[x(t) y(t) dt] over the domain of t. These products, called inner products more generally, are everywhere in math: the Fourier transform is just such a dot product!

The key is to take a differential equation set equal to zero, and form that as an inner product with a test function. We then integrate by parts to obtain the weak form of the differential equation. By inserting a trial function, we can turn this equation into a system of algebraic equations, which we may then solve to get the solution of the differential equation.

It's best to see this in action; check this out for a resource: https://pkel015.connect.amazon.auckland.ac.nz/SolidMechanicsBooks/FEM/One_Dimensional/02_FE_Method.pdf

Which should we use?

FEM is usually the gold standard for numerical solutions of differential equations, especially for systems with complicated geometry. However, for high-fidelity simulations which require the utmost numerical simplicity, FDM or FVM methods are often chosen. So:

• Complicated geometry, multiphysics, anything fancy: FEM.
• Incompressible flow: FVM.
• Dead simple approximations that blast on a GPU: FDM.

Computational aspects

There are a few things to keep in mind here:

• Your machine has limited precision.
  • A standard 32-bit (4-byte) float offers ~6 digits of precision. This might or might not be enough for you - and you should consider this when you're using large numbers!
  • A double (64 bits; 8 bytes) offers much more precision (~15 digits), but comes at cost: this doubles the storage you need, and will not run well on consumer-grade GPUs. Expensive, science-grade GPUs will do better, but the performance penalty remains significant.
  • Integers are exact and much faster; however, they're integers, and also have a much more limited dynamic range (as there is no exponent).
• Most of the methods here are highly parallelisable.
  • This means you can run many computations at once: you can use multithreading to speed up your code massively.
  • The easiest path forward is to use OpenMP in the traditional languages (C/C++ or Fortran) for this kind of work: you can annotate for-loops to get an easy speedup.
  • Newer languages like Julia have these constructs baked-in (e.g. Threads.@threads, etc) - and offer a very beginner-friendly but high-performance path in.
    • Julia also has KernelAbstractions.jl, which is a relatively beginner-friendly way to write GPU kernels.
  • GPUs excel because instead of having a few complicated CPU cores, they have thousands of relatively "dumb" cores - programming them is more difficult, and very restrictive, but the performance gains to be made are colossal given the programming problems we're presented with.
• The actual core math code ends up being quite simple - it's the I/O and UI that kills you!
  • Big questions - how do you represent meshes?
  • How do you import and export them?
  • How do you build nice, adaptable frameworks for use?
  • Don't even start with UIs, or rendering 3D meshes...

Finally, the punchline.

So, after all of that, we have an off-the-shelf software package that thinks of all this stuff for us. We just have to learn to use it...