All right, we're going to start off again with the Lumi Supercomputer and Machine.
Thank you.
Good morning, everyone.
My name is Boris.
I want to talk to you about my adventures,
scaling a finite demand code and the Lumi Supercomputer.
And now we have to improve our machine code for to work on it.
So first thing first, I will describe briefly what is G-Mesh.
Best generator we use in our team.
And then I will describe the kind of simulations we want to run.
And I will focus on two aspects.
The overlap of mesh subdomains.
And now we mix it challenging to have a scaling code on massively parallel machines.
So what is G-Mesh?
G-Mesh is a free and open source machine generator,
which can build 2D 3D meshes suited for finite elements.
And it has scripting language,
binding for mail and widgets,
C++ Python, Fortune, Julia.
It's also a GUI.
You can use it on the common line, so it's very versatile.
And it can generate that kind of meshes.
So millions of elements,
it is multifraid that you can be on complex geometry,
it's very fun tool to use.
And so this is a machinery to use for my own finite element simulations.
My main field of interest is time harmonic wave simulations.
So wave simulations with non-frequency.
It has many applications, such as rig-to-design imaging.
It can be acoustic waves, electromagnetic waves, so it's a broad family.
But each time I want to do finite elements to solve this kind of problem.
And when you do finite elements for time harmonic waves,
you end up with a very large complex linear system with a spasm matrix.
Typically you can go up to hundreds of millions of millions of a known.
At this large scale, the spasm direct software is too expensive to use.
And classical iterative method, such as multi-grids, dumb work well.
So we want to use an iterative software, but this one.
And so this is the reason we use mostly a domain decomposition method.
Once again, it's a very broad family of numerical methods,
which are naturally suited for parallel computing,
because you split your very large problems into many subproblems.
You solve them in parallel.
And you use it as a preconditioner as some kind of iterative procedure
to recover the entire solution.
And so it is an iterative method overall,
but you need to solve many subproblems in parallel many times.
And when you do that, you can solve very large problems,
such as this plane reactor intake,
which requires a billion of a known at least,
because it's a high frequency wave.
So you have many, many unknown points,
but we can make it work.
This particular test case was run on the Lumis supercomputers,
which is a huge PC cluster in Finland,
which has a main advantage of having enough resources
for what kind of problem.
So for instance, back to a plane this time for the exhaust.
It ran on Lumis with 500 nodes,
which is the most you can ask for on the cluster.
And we could solve this problem in 15 minutes,
and solving a problem with 1.3 billion unknown.
In this case, it was 4,000 MPI processes with 16 threads each.
So that was the state of the antinmyteen before I came in.
So it's a high auto mesh split with Metis,
which generates another lapping partition,
and we ran over 4,000 subdomains.
And so that there were two limitations.
It's that going above is a 1,6,000 subdomains.
It's not very effective.
And also it was another lapping partition,
but sometimes you prefer to have overlapping subdomains,
because some solvers rely on it,
and we wanted to know if there were good old bad solvers.
And so these are the two things I want to talk about today.
So for the overlaps, I think a picture is worth a thousand words.
So you open your favorite mesh with GMesh.
You want to split it, say it on four subdomains.
So if I also read the existing, you call Metis.
Provide a partition, and you use it.
We build interface elements.
So you can do the coupling.
Now the new thing is that we can create overlaps.
So we take a subset of the neighboring subdomains,
and we open to them to have an extended subdomain.
And we use this extended subdomain as a computation component
in the preconditioner.
And so this is a new feature we implemented,
which was a very finite venture.
Once it is done, the API is pretty simple,
and you can query the overlaps, all sort of the boundaries.
So the implementation is that we have the view on existing elements
for the overlaps.
But we need to create new elements for the boundaries.
In particular, we have two kind of boundaries.
We have a outer boundary, and what we call inner boundary,
which may be not completely appropriate.
So the green part here is the extension of the already existing boundary.
And the red part is the artificial boundary of the overlaps.
But it actually is a curve inside your computation domain.
And so it's critical to know the difference,
because you want to apply different boundary conditions on these two families of curves,
because it impacts correctness and convergence of your thing.
And so we added these features into GASH.
And once it's done, we can do some numerical experiments
to perform large scale simulations on the main.
This will be my test case until we end the talk.
It's a geophysics test case for aquistics simulations.
So it's an heterogeneous medium with varying wave velocity,
blue with water, where the waves are slow,
and yellow red are rocks or soil,
where waves travel much faster.
And so it's an interesting problem,
because you have varying wavelengths, reflection, reflection,
and interesting physics happening.
And I want to simulate ways at higher frequency possible.
And higher frequency, the more expensive a problem,
higher frequency means shorter wave lengths,
which means a finer mesh.
And so more degrees of freedom.
And so I go to from 1.6 to 6.25 hertz,
which for in our mesh translates from 18 to 600,
11 to 600 million of a node.
And we do it in a weeks scaling way.
So we want to keep the problems ice constant,
and we add more resources as we increase the frequency.
And so we need one process for every 85,000 unknown.
And so using the resources from Lumi,
it went from 2 to 128 nodes,
each node having 100, 28 CPUs.
I mean, threats.
And so typical example of a kind of mesh we have.
We have a lap is not shown, but it was computed.
So you can see each color is a partition,
so a different subdomain.
We are designed to have the same size in terms of element count,
followed balancing, but we actually represent different sizes
in terms of geometrical space.
And because we have smaller elements in some areas,
and larger elements in some areas,
to match local wavelengths.
In the original 5 format we used,
we have one file, passive domain.
Each file contains the necessary elements, the nodes,
so this part of the mesh.
But we chose to duplicate the mesh topology.
So each file contains information about which subdomain,
to choose which subdomain.
Once the name of interface, between the subdomains,
and that kind of things.
Because it was simpler.
But actually at some point it will turn out to be an issue.
I'll come back to it later.
But now that we have overlapping subdomains,
you also need to export elements from the overlaps.
So the files get a bit bigger.
You also need more bookkeeping information
to know which subset to use and what kind of things.
You also need to store boundaries of overlaps.
So the whole thing is a bit more expensive.
And so we start with a fairly naive implementation
of exporting all elements from our subdomain,
and the elements from the neighboring subdomain.
Which is a bit too much because we don't need all elements
from the neighboring subdomain.
But it's simpler and it's just a mesh.
So we don't expect it to be too expensive.
So there is some redundancy.
And it worked on small problems.
We managed to run up to 40 million unknowns.
And then going above, it's getting complicated.
Because we have too much duplication.
And sometimes you can load all the meshes
because you have way too much data.
And especially as you increase the problem size,
we have a more regular mesh,
which translates to more of a laptop domains,
so more elements to export.
So we need to do something a bit more efficient.
And so we decided to export only the necessary elements
in the necessary nodes.
Which takes some time because you need to loop
of elements to see if you are needed or not to keep a set.
But it makes for smaller files and mesmemory consumption in less ways.
So we do that for the elements and for the node.
And it helps.
We could go a bit higher and solve 80 million non-problems.
But going above, still crashing.
Because we run out of memory.
And so as someone who loves to waste energy,
just for a render job with more memory to see if it worked.
And it didn't.
And I was a bit confused because we export only what we need
or subset of elements, subset of node.
And it turned out that the mesh topology I mentioned.
But we keep data on which the domain is connected to.
What's the domain and what kind of things?
Which is deprecated across every file.
When you have a big thousand partitions, it's not negligible anymore.
And that actually makes up for most of a memory you consume.
And so what we ignored because of the minor overhead.
It was actually wasting more than half of our memory.
Just to use this information.
And so we felt a bit dumb when we noticed.
And it was a bit annoying to fix because many parts of the code
relied on this assumption.
So it took almost a month to fix.
But once we did it and everything worked well,
we could go much higher and we could solve problems with.
200, 300, 1600 unknowns and go to up to 6000 subdomains.
Which is four times what we could do before.
And it's with other lab.
And we could even go I go in theory.
And now the mesh fans are much, much smaller.
And so we reached the efficiency we were hoping for.
And using this, we could publish paper, comparing other labping and non developing numerical methods.
And it turned down to non developing one is a bit faster and requires less memory.
So all this for that.
But it was worth investigating.
Very small subset of our paper where we do some weak scaling studies.
It's a bit complex to read if you're not familiar with a problem.
Because we expect that even in the case of perfect weak scaling,
the timing increases with the frequency.
Because the cost per iteration is stable.
But the number of iterations still grows.
So we expected the linear relationship between time and frequency.
If we had resources as a frequency increases.
And we could show the oras on the left is overlapping method.
And the oras on the right is non developing method.
And non developing method is faster, shorter and longer.
So once we publish this paper, we could publish a bit of implementation.
So we recently made the overlap boundaries,
but with high order elements.
Which changed the way we build the overlaps.
Because it was a bit naive on the first time.
Instead of having a subdomain looking for elements on its neighbors.
Each subdomain looks at its boundary and tells the of a subdomain.
This should be a overlap.
And it's a much simpler and much faster.
And we finally merged.
Last week all these changes in the main branch of gmesh.
It's not that part of the release, but it will be in the next one.
And so if you're interested, you can already try it out.
So it was a bit story of a big story of attack on in my life last year.
It was a long adventure.
So it's a good lesson on optimizing and scaling things.
Because at every step, a new issue arrived.
And it's never the one you expect.
So it was a lesson in always profiling, understanding assumptions.
Because sometimes you think you know when you have things consuming too much memory.
And it's actually something completely different that you didn't think of.
But it was very interesting to always improve iterate, eliminate waste at each step.
And going further forces you to find out what's the main issue.
And it's very fascinating things to to live.
The last thing we haven't fixed yet is that we can do fast simulations when you load the mesh.
But the time we take to export the mesh properly is a bit long.
Because we do extra checks to ensure we only export what we need, everything is correct.
So it could be optimized a bit further.
But we came a long way since last year.
That's about it for me.
Just a small advertisement that we were seeing the first gmesh user meeting in the edge of the summer.
So if you happen to be a gmesh user, feel free to come by and discuss this kind of things.
And I'm open to questions if you have some.
So the question is why I didn't use an adaptive mesh measurement, right?
Actually we did.
No, it's not adaptive.
It's a, um, to say it.
We create a mesh in advance.
But we know where it has to be fine and why it has to be because because it's from the local wave speed which we is known.
So once you have partitioned your mesh changing the mesh density would be a bit of a mess and not sure we plan for it.
But in this case we already know the ideal refinement level.
Okay.
So the question is does gmesh support GPUs?
Um, great question.
Working progress.
So the machine process is itself still wants on CPUs because it's not very GPU friendly.
But one of my colleagues is performing the same kind of simulation I'm doing.
But on GPUs.
So it's not very natural because it's pass on GPU everywhere.
And it's not naturally suited for GPUs.
But essentially you do something as I do with subdomains but you do much smaller subdomains.
And you instead of having one subdomain per thread.
You have maybe 100 and 1000 subdomains per GPU.
And you can use a 10 or 100 GPUs if you want.
Which is one more reason to partition efficiently.
And in the future we will need to have maybe 100 subdomains per file instead of one file,
because each PC and means don't like when we have 16,000 files as input.
So far we haven't murdered me yet.
So the question is, which GPU we used?
Um, which technology?
And I have to ask my colleague for details.
So we have access to a Belgian cluster which is Nvidia GPUs.
And we have access to Lumie which uses AMD.
We mostly prototyping of a smaller one which is Nvidia.
But we still work in progress.
But it's we're trying to be compatible with both.
But there's still some compatibility issues with some libraries.
And we're comparing different backends for the algebra part.
Because we have some dense backends and sparse backends.
And it's very exploratory.
So the question is, did I consider out of core competitions?
Do you mean for the simulation part of the...
Yeah, so for mesh part, it's not really an issue for now.
The issue is time to export all the files.
But we can do it.
It did a lot of memories.
Sometimes you use what we call a fat node and the cluster.
So we just more RAM and less CPUs because we just need a lot of RAM.
So maybe it can go up to a terabyte of RAM.
And for the competition part.
The solvers we use aim to avoid that.
So you can use sparse direct solvers, such as memes which have an out of core capability.
But it works with its slower.
And even if you have enough memory, it's actually like five times slower than what we do.
And so we save memory with the time.
And we don't have to mess up the disk.
But it depends on the application.
Because we're investigating case.
We solve for multiple right-hand sites.
I'm focused on case with 50 or 100 right-hand sites.
But if you have thousands of right-hand sites,
maybe it's worth it to have a factorization of your matrix.
Yes.
So the question is what is a gmesh format like?
So it's a custom format, but it has a binary and text version.
Binary is just a bit faster to read.
But we all have the implementation support both ways.
And when you open the file, it detects which one was used.
Can you repeat the last part?
Okay, so the question was,
where if we use GPU partitioners,
and come again?
Okay, yes.
So the question is about alternatives to matrix?
For now, we only support matrix.
We thought about parameters,
but slightly sensing issue,
and the gmesh itself is still controlled.
So it's not immediate.
We are working on prototype times of all main partitioners,
based on space feeling curves.
Which has supposed to be much faster,
but if very rough partitions.
Not very familiar with GPU ones,
but it would be nice at some point probably.
And we've regard to the physics,
there's no particular trick.
We just give a graph of volume element,
so it may just know the connection between T-Trader,
and we just ask for partitions with a similar bit of T-Trader
in each chip domain.
Last one.
Exiderer?
Yeah.
So the question was,
is there support for exiderer or over candle elements?
Very support for exiderer elements.
Very recently, like last week,
they started working on prototypes for political elements.
But it is very much of a draft for now.
But exiderer, yes.
At least there exist in the format.
Not sure we algorithm for generating them as advanced,
but you can exist.
Alright, let's go.