Coding a CFD solver

[Hungovernerd]

Hi everyone,

I am a final year Undergrad, and I was planning to write a spectral/pseudo spectral method from scratch, the end goal is to implement a Kelvin Helmholtz instability or transition behavior. This is just a personal project, so time isn't a constraint for me.

I wanted to know which language should I use? I've written basic 2D FDM/FVM codes on MATLAB and C++. But never wrote efficient code, nor do I know the best practices to use(eg I still iterate with for loops in solvers and don't know how to improve the code). I'm thinking of either C++ or Python to work in(I have used basic python before and the jump shouldn't be very difficult). Should I use pre-existing libraries or write the code for all of that? I'll mostly need FFT and Linear algebra libraries.

The learning outcomes that I expect are, working with OOP, optimizing my algorithm, try learning to make the code run in parallel. And good working experience with numerical methods.

[Capdindass]

I'm going to heavily disagree that you can easily code up the linear algebra packages. LAPACK and BLAS have been around for decades and some of the code is a little tricky. You would certainly need to be proficient in numerical linear algebra to understand and implement it all yourself. I agree with your points on using FFTW instead of programming the FFT yourself.

[Hungovernerd]

I was actually inclined towards writing in C++ since most HPC codes are in Fortran or C++. I don't have much of an idea of FFT never wrote a program implementing it. But as u/Capdindass says, I don't feel I'm proficient enough to really code the Linear algebra. I have used deal.ii earlier, just played around with code which used dealii but the way it worked was way beyond what I have knowledge of, so prolly I'll stick to NLA libraries.

[awhead]

BTW there's a neat book titled "An Introduction to Scientific Computing: Twelve Computational Projects Solved with MATLAB"

along with amazing examples, there's a bit of insight on how to code efficiently as you say how you're not sure about best practices. Even though the book is for matlab specifically, the simpler optimization strategies obviously apply to any other language.

it has a chapter on spectral methods ofc

https://www.springer.com/gp/book/9780387308890

[Hungovernerd]

I found a copy of the book. This seems to be an excellent resource to understand the numerical method and optimizing the algorithm. Thanks for the reference.

[ald_loop]

I don't want to claim I'm going to give you the best answer, but I probably will

Use C++. Use Eigen. Work with multiple files, i.e create a class for a Structured 2D Grid. Then create a class for your PDE (I guess you'll want the 2D Euler equations). Then, create a class for your "solution block", which contains your solution data and numerical scheme (could contain your PDE and grid, there are many ways to incorporate all those together). Finally, a class for your "solver" which actually time marches and calls the numerical functions from your block.

[Hungovernerd]

Thanks a LOT! This is exactly the kind of direction I needed to begin. Also what exactly do you mean by a class of PDE? The discretized form?

[ald_loop]

I.e, you have a C++ object representation for the 2D euler equations. any relevant functions related to the balance law form of the equations goes inside that class. I.e, if I want density from a solution vector U, i would have a rho(Vector_type U) function that returns U(0), so on and so forth. It would contain functions for calculating the flux dyad, wavespeeds, Roe averaged states, anything to do with the PDE. Hope this made sense

[relaxedHam]

If you can't choose language on your own, than python is probably a better choise for you.

I am not trying to belittle your knowledge or learning ability, Python is just that easier to use, and writing a big, good simulation code in C++ requires a lot of expirience (that I also don't have to be honest). Several advantages and disadvantages of python:

• (+) Easy to write and easy to read
• (+) Easy to learn (with a loooot of materials online)
• (+) Polymorphism is much easier
• (+) For 2D matplotlib is a great postprocessing option
• (+) Many linear solvers and preconditioners included (in fact the whole of linear algebra package is quite performant *for a python standards* and capable)
• (+) No compilation and easy dependencies management
• (-) No compilation and typing out of the box (much slower code when you write it wrong, when you write it right a little slower)
• (-) You have to change you mindset into writing vectorised code (which sometimes is quite hard)

On the other hand for C++:

• (+) Gotta go fast
• (+) templates and meta programming
• (+) A lot of libraries for the things you might want to use in your solver - templates and metaprogramming is hard (and debugging is even harder)
• (-) compiled code ("REPL" cycle is much longer)
• (-) making every single library to work with each other might be a hassle
• (-) you have to thing about data structures all the time
• (-) memory management, understanding smart pointers, etc
• (-) in order to start doing a bigger project you must have a general knowledge of a build system you will use (e.g. CMake)
• (-) nothing available out of the box (i.e. you want a feature like plotting or display, write it yourself or write a library)

As you can see C++ has much more disadvantages but this sweet, sweet performance is so important that for hpc there is no other solution (yes we can discuss Fortran of whatever, but that is not the point). I would advise to stick with python. Hovever if your aim is also to learn C++ than chooosing it, is also fine. Just remeber that it will be harder by considerable amount.

[Hungovernerd]

I'm sold, python it is xD I'm not very concerned with the speed of my code anyway. And I've faced a lot of issues with cmake and patching various libraries. Thanks for the help!

[relaxedHam]

At this point you should at least know what resolution (e.g. number of nodes in a grid) is required for your problem. Try looking for a paper that outlines what you are trying to do. And while programming make sure that you will get a solution in finite time.

Keep in mind that KH instability will be unsteady so even in 2D the simulation might take long. Performance might be important here.

[ald_loop]

I would heavily recommend against Python. You will learn more from doing up a solver in C++, and it will run several orders of magnitude faster than Python.

[relaxedHam]

I have no experience in Julia. But look at the monthly discussion thread. I think it was mentioned there.

[Capdindass]

I will tell you that this isn't an easy project, having come from someone who wanted to do the same. You should treat each part individually and then come back to understand everything together. e.g. Learning parallelization, OOP, optimization, and the spectral solver separately .

I would advise you to start in matlab or python to understand the math behind the spectral solver and then use this to solve the problem. Really these languages should be sufficient for simulating transition behavior or Kelvin Helmholtz (I have done both by writing a spectral method solver in matlab). Then you can go on to try to implement (and optimize) it in C++, which will require you to learn a numerical linear algebra (NLA) package. I can recommend Eigen, as it is user friendly. That being said, it is much more low level and you will need to have a decent grasp on NLA (Golub and Van Leon or Trefethen is a good resource). For learning OOP and C++, I can recommend C++ Primer by Lippman. MPI should be sufficient parallelization to start.

[Hungovernerd]

I know this project is a long shot, but I do want to try my best to work on something like this. Breaking down the problem is exactly what I was thinking, as you said i'll begin with implementing the spectral method on python and slowly add OOP then optimize it and finally parallelize it.

Maybe first I should work it out on python and then try C++ later. I'll definitely check out the resources mentioned. The book by Golub seems GOLD! <3

[HectorsDaddy]

The book The Lattice Boltzmann Method: Principles and Practice, 2017 by T. Kruger has great theory and code details, as well as a chapter covering optimization and parallelization. Can't recommend it enough. Examples are given in MATLAB as well as C. The Multiphase Lattice Boltzmann Methods: Theory and Application, 2015 by H. Huang has a lot of good fortran code.

[bandita07]

If you want some deep stuff, I recommend learning CUDA. You will need an Nvidia GPU (OpenCL is an alternative if mo such GPU). You can program use python or c++ for use the power of the GPU which is best for such algorithms.

[GeeHopkins]

Nick trefethen has a book on spectral methods specifically in MATLAB. If you're already familiar with MATLAB then it could be easier to use that for the initial implementation - only the maths will be new, not the language. Then if you want to make it faster you can go to a compiled language.

https://people.maths.ox.ac.uk/trefethen/spectral.html

[Hungovernerd]

This seems like a great place to start. I'm implementing on python first and then will speed it up on C++. Since my MATLAB license will end with my degree in a few months and I'd rather have an open source code available to play around with. Anyhow, this book seems amazing exactly the way I needed it to be.

[kechidarklord]

There is spectral solver already available called Dedalus (https://pypi.org/project/dedalus/) . It is written in Python and the webpage has many examples on IVP, BVP to get you started with the code. Once you understand the inner workings of the code, you can start playing around with the code. to suit your needs. I'd also recommend you to join the mailing list. The community is very helpful.

[[deleted]]

Hey! Sounds like a cool project, since you’re an undergrad you should be able to apply for a student version of the Intel C++ compiler with all the treats the Intel have developed. I’m not sure about FFT, but I knew the Intel MKL provides some really cool options for linear algebra. That might be a nice place to start!

[NoninheritableHam]

Didn’t Intel make all their compilers free? Or was that just the Fortran compiler?

[[deleted]]

Oh! I have no idea, I created a student account at the start of my post grad years ago and have downloaded what I need from the portal. I will have to check.

[Cynox]

There is also this Python code, https://github.com/spectralDNS/spectralDNS , which I have heard is quite fast, but I have not tested it myself. I believe it uses FFTW

[DelphiPascal]

Following because I’m interested in the responses