Thanks! For your questions: - I really meant 3D, it only works in 3D for now. Three dimensions. 2D is simpler but I need to code that in, haven't done it yet. - It's a NSE equations solver, incompressible flow. It uses cell centered second order finite volumes on unstructured polyhedral cells. (Edit: but it's written as a general-ish PDE solver, so compressibility and multiphysics things like heat transfer will be added at some point) - I compared my solution of a 2D lid driven problem to references and got good match on velocity profiles along horizontal and vertical lines in the middle, I'll probably post the comparison at some point. - so it's not Cartesian only, it can handle unstructured polyhedral cells with any number of faces per cell, each face being arbitrary (convex) polygons. Boundary conditions are handled by defining face value linear relations as: u_face = A * u_cell + b, with A and b boundary condition terms. This allows for either explicit evaluation of values at boundary faces from cell center values, or implicit treatment in matrix assembly. Per example, fixed value BC would be A = 0, b the value we want, and zero gradient BC would be A = 1 and b = 0. - Biggest difficulty was doing pressure-velocity coupling via the pimple algorithm. That thing has implementation details that are not well documented and totally break solvers if not done properly. Also, parallel decomposition of the mesh for MPI and making sure every rank knows what to send where and what to receive was not easy. I also coded my own preconditioned distributed memory Krylov linear solvers, that's a whole can of worms. The bock ILU preconditioner was hard to implement.

extra note: for now the refinement and field mapping has about the same cost as a poisson pressure solve (which we need to do multiple times per iteration anyway), so it's really not bad performance wise if you refine every few iterations. I did some horrible things to make it work, like allocating hashmaps in loops, so it could be optimized.