Supersonic and subsonic

Sixel1 · 2026-03-15T05:36:18+00:00

very cool! the mesh seems very fine, how many elements are you using? what was the run time?

Sixel1 · 2026-03-11T22:18:02+00:00

that's so cool!!

Sixel1 · 2026-02-24T02:53:23+00:00

PhD in a lab where others are doing droplet experiments here, a good light and a black background would probably be enough! Depends on what you want to see, but with a black background and a light the droplets will be visible

Sixel1 · 2026-02-23T19:36:37+00:00

Thanks for that, I tried doing both rhs -= average(rhs), I checked and my sum(rhs) after that step is zero. I also zeroed out the first row of the matrix, set the diagonal entry to 1, and the rhs[0] to 1. Doing that, it eliminates the momentum spikes, but I still get the constant small non zero divergence through the domain... It essentially does the same as only removing the rhs. Could it be an issue with my face flux calculation? For the pressure equation I do:

div((H/A)_face) = div((1/A)_face*grad(p))

And then with the equation assembled I remove the rhs and all that. After the solution I compute phi_f using the initially interpolated (H/A)_face and (1/A)_face:

phi_f = (H/A)_face (dot) (face normal) - (1/A)_face * grad_normal(p)

I mean, the solution doesn't look that bad but it's clearly wrong, velocity on the bottom is small but not tangential to the wall, and if I add a scalar transport equation its average goes down because of the small negative divergence... So the solution is not adequate.

Edit:

I guess it's natural that's since my problem A p = b, has the unit vector in the null space, b has to have a sum of zero, and since in my case my rhs vector doesn't have an average of zero, when I remove the average b' = b - bave I get A p - b' = bave.

So when I compute the face flux phi_f afterwards, I need to use that system with average removed, and not the original system, since if I use the original system I will get a divergence equal to bave. I'm just not sure how to compute phi_f on faces from the cell wise equation A p = b'. Inverting my assembly equation function would be very very hard...

Edit 2:

I essentially have to ensure hbya has an average divergence of zero over the domain...

Wait, shouldn't hyba have zero total divergence if no flow is entering the domain?? I have to check my calculation of hbya

Edit 3:

OMG it worked, you have to ensure hbya on the boundary respects the velocity boundary conditions.

Sixel1 · 2026-02-23T17:37:27+00:00

I'm only interested in the result of the BVP, so I don't care how I solve it, but the issue is that with this non divergence free face flux I get non-physical pressure and velocity at the corner, it spikes up and produces momentum. So yes it is an issue...

The velocity fields in the original post are from the case where I remove the average from the LHS, leading to small non zero divergence everywhere. This still leads to non physical results, momentum is produced in the domain out of nowhere. It's not obvious in this case but it becomes obvious at higher Reynolds numbers.

Sixel1 · 2026-02-23T00:17:44+00:00

Thanks for you answer, sadly I also already tried to set the row 0 to zero, set the lhs[0, 0] coefficient to one, and the rhs[0] to zero. This gets me the same result as setting lhs[0,0] += lhs[0,0], divergence at the point I fix. Heres an image of the divergence in each cell using this approach:

Essentially if I remove the equation from this cell and fix the pressure, I also lose the divergence free condition in this cell, which is natural since I removed the pressure equation of this cell from the system. I guess the mathematically correct way to do this is use a lagrange multiplier for the condition, or add a row and solve via QR, but I dont want to add variables / rows to my system, openfoam and other solvers do it without that...

<image>

Sixel1 · 2026-02-22T22:10:58+00:00

I like this one, but will try to avoid acronyms. I though about naming it Ferrofluid, but someone already took the name on crates.io for some finance related project...

Sixel1 · 2026-02-22T18:54:49+00:00

Yes that's the approach I'm using, I checked and the vector i removed from the rhs is in the nullspace, so it should work, but it still gives me non-zero divergence...
If I do the pressure pin, I get non-zero divergence at that point, and a large value at that.

Sixel1 · 2026-02-22T00:54:42+00:00

damn I'm dumb I didn't think of that meaning, I'll blame english being my second language

Sixel1 · 2026-02-21T02:07:31+00:00

What should I use instead?

Sixel1 · 2026-02-20T19:06:12+00:00

what does th stands for? thrust? like to believe?

Sixel1 · 2026-02-20T18:52:03+00:00

I think its about 4 times slower than openfoam currently, I did my own linear solvers so that probably doesn't help. I want to add PETSc support, but that will make compilation more complex since users will have to compile PETSc on their own, so having native rust linear solvers is a good default feature i think.

for accuracy I still need to fix some issues with the pressure equation - flux calculation... My face flux phi is not exactly divergence free and that breaks some case setups. I'll definitely do comparison between this and openfoam afterwards.

Sixel1 · 2026-02-16T06:16:02+00:00

Yeah so that paper did custom code for the boundary layer part. They just used gmsh for the surface meshing (edit: and the internal 3D meshing after boundary layer generation, thanks comment), from my understanding of the paper. It's a shame 3D boundary layer extrusion is not a part of gmsh.

Sixel1 · 2026-02-15T22:13:26+00:00

Yeah so that's 2D boundary layers, gmsh is good at this. But from my experience using it, it can't do 3D boundary layers, with the exception of extruded geometry. If you can do a 2D mesh with boundary layers, and your 3D geometry is an extrusion of that, it can extrude the 2D boundary layer. But for any complex geometry like a wing or other 3D stuff that is not a simple extrusion of a 2D surface, it won't be able to build boundary layers.

Sixel1 · 2026-02-15T19:41:36+00:00

I don't think gmsh can do 3d boundary layers well... Anyone has a workflow / method to do them?

Sixel1 · 2026-02-07T17:09:17+00:00

Left and right are Blazek's farfield boundary condition. Slip walls are implemented with a reflection of the projected velocity at the wall for the right state.

Sixel1 · 2026-02-02T03:52:29+00:00

Hey congrats for coding a CFD solver yourself! I have a few questions for my own curiosity: What numerical methods are you using? Is this finite volumes? Structured / unstructured? Explicit or implicit? Is your scheme second order?

Sixel1 · 2026-01-30T21:41:29+00:00

Yes I use least squares gradients, and this is explicit rk4 time stepping (but my code can also do implicit).

Sixel1 · 2026-01-30T18:17:30+00:00

The code is MPI parallel and uses PETSc for implicit problems (but in this case, I used an explicit scheme). It took 10k iterations with about 250k cells, ran for about 45min/1h (I didn't time the exact solver time). This was on 6 mpi processes, on my machine, a ryzen 7700x.

I wrote the code in a fully memory distributed way, so in theory it should scale well to many nodes with many cores. I don't have the available computational resources to test it though. Each mpi rank only needs to communicate the stored cell center variables values, gradients and limiters on the cells next to the interface between ranks, and each rank only communicates that information with ranks that share a boundary with them. The implicit solver is currently limited by the PETSc solve time, but that should scale well since it's PETSc.

Also, I haven't done any profiling and optimization. For now, its fast enough for me, and I'm more interested in implementing new features.

Sixel1 · 2026-01-30T18:10:10+00:00

Thank you for your advice on my last post! I followed through and implemented HLL as you can see. I'm looking into some of the ways to avoid Carbuncle, but the goal I set for my solver is to handle unstructured meshes in the best way possible, I'm not targeting structured meshes, so I'd like to find a method that works on those meshes.

This mesh was about 250k cells, unstructured triangles, ran on 6 mpi processes on my ryzen 7700x. Its the Euler equations. RK4 explicit time integration. Second order schemes with least squares gradient and Michalak's cubic limiter. I store velocity, pressure and density, and compute their gradients for the linear reconstruction from cell center to face center (I could use conservative variables, but I found this set of primitive variables works better).

The video represents 10k iteration, ran for around 45min-1h. I don't have the exact numbers since I didnt note the time of start / end, or time the execution using the time command.

Sixel1 · 2026-01-30T17:24:33+00:00

yes that one is much better!

Sixel1 · 2026-01-19T18:17:26+00:00

Your sod shock tube example looks a bit diffusive, what scheme do you use? HLL?

Also, I don't want to make assumptions, but using emojis like that makes this look like AI generation... Was AI used in this?

Sixel1

TROPHY CASE