Periodic FE space in parallel computation (PETSc)

Dear FreeFem developers and users,

I’m writing since I’m facing a trouble conncected with the definition of periodic finite element spaces in case of parallel computations.
I am just trying to define a very simple periodic fespace in a cube and everything works properly (or at least without any error) if I define it in a sequential problem (or even using parallel code structure, but using a single processor), while if I increase the number of processors adopted in the simulation I get the following error:

Exec error : Periodic 3d: the both number of face is not the same

I have the suspect that this is due to the mesh partitioning which is introducing “new labels” given by the “surfaces” along which the mesh is fragmented, but I’m not able to guess how to solve the problem!

Is someone that have already faced this issue? How did you solved it? Is it possible to solve it?

I am already thanking a lot whoever will help me!

sincerely yours,

Marco

P.S. just for completness I’m attaching here the few lines which I am using for the definition of the above mentioned finite element space in which I find out the problem!

load "msh3"
include "macro_ddm.idp"

mesh3 Th = cube(10,10,10);
if (mpirank == 0)
{
  cout << "Mesh labels before partitioning = " << labels(Th) << endl;
}

load "PETSc"
macro dimension()3// EOM            // 2D or 3D

macro defPart(u)u// EOM             // partition of unity definition
macro initPart(u)u// EOM            // partition of unity initialization

int[int] n2o;
macro ThN2O()n2o// EOM

buildDmesh(Th) //now the mesh is local

cout << "Mesh labels = " << labels(Th) << endl;

func Pk = Edge03d;
fespace vh(Th,Pk, periodic = [[2,y,z],[4,y,z]]);

See FreeFem-sources/diffusion-periodic-2d-PETSc.edp at master · FreeFem/FreeFem-sources · GitHub for how to use periodic boundary conditions.

Thanks a lot !
This is exactly what I was looking for, but I was unable to find it out!

By the way I would like to ask you another question.
Imagine to use this strategy in the case in which the domain is a spherical wedge with opening angle which is relatively small.
If you impose as periodic surfaces the two planes forming the spherical wedge what you get as output is a mesh decomposition in which there is one of the “submeshes” which includes a much larger amount of cells with respect to the others since it must include both the planes used as “periodic surfaces”.

Having not an even splitting of the number of cells between the different processors is not efficient from the computational point of view.

So does an optimization strategy exists to improve the mesh decomposition or is it better to enlarge, if possible, the angle of the spherical wedge trying to obtain a more eaven distribution of cell number between the different processors?

I hope that my question is clear!

Thanks a lot!!!

Your question is very clear. Partitioning meshes for periodic finite elements is non-trivial. Often, the best choice is to come up with an ad-hoc solution depending on your geometry and the problem at hand. One such solution is exposed in the other example FreeFem-sources/diffusion-periodic-balanced-2d-PETSc.edp at master · FreeFem/FreeFem-sources · GitHub. Of course, you’ll need to adapt this to your problem. I can help you with that if it’s not clear.

I try to study the example you are mentioning trying to adapt it to my geometry.
I will contact you again in case I will not be able to understand it on my own!

I really want to thank you again for the precious help!

Dear @prj,
I have spent quite a lot of time trying to understad the solution that you proposed to me for the optimization of the mesh partitioning in case of periodic finite element spaces.

The first step which I have followed is apllying the strategy you suggested to me to a different geometry (a very simple square) just to understand if I was able to get the point. Indeed I have wrote the first script I am attaching (PeriodicBC_simple.edp).
PeriodicBC_simple.edp (5.6 KB)

On this topic I have a question:
Making use of this strategy if you increase the number of processers you are using you reach at a certain point a condition in which one (or more) of the mesh part is not in contact with the periodic border. This is not good for the strategy adopted since you get a bug (you can check it simply running the script I have attached with 8 or more processors). Do you have solutions to this kind of bug?

After that I have decided to check if the partitioning I was able to develop was suitable to the solution of a simple problem (I have proposed a simple thermal diffusion problem). I have developed the sequential benchmark in order to check the solution (testTermicoSequential.edp)
testTermicoSequential.edp (2.0 KB)
Getting out the reference result:

I have subsequently moved to the solution of the very same problem, but using the partition I have developed above. At this point I have founded out a couple of issues:

  • The results can apparently be obtained only using a certain number of processors. Indeed using, for example, 5 processors the reference result can be reproduced:
    , and this is true also using 3 processors. By the way changing the number of processors to 4, 6 or 7 the PETSc error “Argument out of range” appears. On the contrary using 2 processors the result can be obtained, but the distribution appears to be disturbed
  • Another issue has been founded out using the restrict command to “move” the results in the non fragmented mesh (line 211 of testTermicoParallel.edp). Indeed the restric command causes a segmentation violation which I’m not able to understand.

Here below the parallel script I have used:

testTermicoParallel.edp (7.4 KB)

That’s what I get out trying to implement the strategy you suggested me. I would be extremely grateful if you’ll help me to understand what I’m missing.

I thank you a lot already!

Yours,

Marco

Do you have solutions to this kind of bug?

I don’t have any bug with your PeriodicBC_simple.edp with 8 subdomains.

With respect to testTermicoParallel.edp, there is no support currently for the use of n2o + buildWithPartitioning, so the segmentation fault you get is the expected result (you can check that the n2o array is empty). If you know a priori the shape of your domain, you could define the partitioning near the periodic surfaces, which would make the partitioning more robust even with large number of processes. The problem right now is that as you increase the number of subdomains, the partitioning gets worse (badly-shaped subdomains), and then we fail to build the distributed fespace.

First of all thanks again for the prompt answer!

I really do not know which could be the reason for which you are not finding out the bug in “PeriodicBC_simple.edp”. What I get is what follows:

image
and this happens at the point in which the script reaches the line 45 for ii = 4 since it is trying to truncate a mesh following a boolean expression which is always false due to the fact that no one of the cells in the considered subdomain is close enough to the domain boundaries (this is easier to understand looking at the following picture)


Isn’t that happening in your try?

For what concerns the use of n2o+buildWithPartitioning, do you have any suggestion to get the same results, and so having this “movment” of the results on the “backup” mesh? This would be useful to me for postprocessing purposes.

To conclude, I know the shape of my domain, so can you suggest me some example in which I can learn how to define a partitioning in an arbitrary way and not in the “automatic” way?

Thanks again!

Isn’t that happening in your try?

No, I have a slightly different mesh partitioning for n = 8 and it seems that all subdomains are touching the boundary. But you are right that ones has to be careful, it is illegal in FreeFEM to do Th = trunc(Th, false), so you have to make sure that the Boolean expression in your trunc is not uniformly false.

Instead of n2o+buildWithPartitioning, I would use n2o+UserPartitioning, cf. https://github.com/FreeFem/FreeFem-sources/blob/master/examples/hpddm/diffusion-cartesian-2d-PETSc.edp#L28. I think it’s already working, but if it’s not, I’ll add it (on the other hand buildWithPartitioning is an old interface and I have no intension of updating it).

I know the shape of my domain

What is it? Always a square? In that case, I would just extract the boundaries and one or two elements that are adjacent to it, like mesh ThNearBorder = trunc(Th, x > 0.9 || x < 0.1 || y > 0.9 || y < 0.1); this mesh would be assigned to one or two processes, depending on the size of your square, and then, for the rest of the mesh, you can do a “normal” partitioning using METIS with partitionerSeq.

Ok! That’s fine. I’ll try to find out a way to avoid having uniformly false boolean expressions in the trunc command!

Thanks for sharing this example, I will immediately check it!

For what concerns the domain I am using a square only now to understand how this procedure works. Once I have understood how to make it works,I’ll need to upgrade it to a 3D domain. Indeed I’ll finally work on a sphere wedge.
In any case how is it possible to assign the ThNearBorder to some processes? How to avoid that other subdomains obtained with partitionerSeq would be assigned to the same processors as above?

See this other example: FreeFem-sources/restriction-2d-PETSc.edp at master · FreeFem/FreeFem-sources · GitHub. In your case, you don’t want to call partitionerSeq on both part of the domain, just on the interior. For the elements near the periodic boundary, that’s where you want to do potentially some decomposition by hand.

1 Like

Here I am again @prj … I really apologize for disturbing you a lot, but probably without your help would be impossible for me finding out the proper, and hopefully, definitive solution to my problem.

For what concerns this suggestion:

I succeedeed in completing the partitioning with the new standard (UserPartitioning), both using a “part” generated by partitionerSeq and using a “part” which I produced taking inspiration on what you did in the above mentioned example too.

by the way the execution is freezed at the point in which I add the “createMat” command without reporting any error (even with very high verbosity). This behaviour is very hard to be understood to me.
I am attaching here below the two scripts which are generating the situation I have presented above.
testTermicoParallel_v2.edp (6.9 KB)
testTermicoParallel_v3.edp (6.2 KB)

do you have any idea of what could it cause this “bug”?

I really what to apologise for the extremely long thread and thanking you a lot again.

I don’t really understand what you are trying to do (with the latest approach I suggested, you are not supposed to have to deal with movemesh and stuff, just the initial domain that you split in two and that your partition with partitionerSeq on one hand and yourself on the other hand), but in your script, the problem comes from the fact that you modify Th with a trunc after buildDmesh and then calls createMat. That’s not supposed to work, you should not modify the input mesh of buildDmesh yourself.

I’ve adapted restriction-2d-PETSc.edp to show you what I meant originally. See attached (must be run with more than 3 processes).
restriction-2d-PETSc.edp (2.8 KB)

The movemesh practice I have included is linked only with the fact that, actually I was thinking it was needed to satisfy the periodic fespece I am using. Is the solution you are proposing in restriction-2d2PETSc.edp working on its own in case of periodic fespace? Shouldn’t the periodic sides communicate in some way one with the other (this was actually the reason of the movemesh part if I have understood properly)?

Anyway… tomorrow morning I’ll try to implement this “domain splitting” that you are suggesting and I’ll let you know!

Thanks again!

Shouldn’t the periodic sides communicate in some way one with the other

Yes, of course. In the previous .edp, periodicity should only be applied on the top and bottom sides. If you need periodicity on all sides, then something like this will do the trick:

partCondFE = ((x < 0.25 || x > 0.75) && (y < 0.25 || y > 0.75)) ? 0 : 1;

Now I have understood! I succeedeed in increasing the number of part in the “close to boundary” region too!

I’ll try now to pass to 3D hoping not to disturb you anymore!

Thanks a lot @prj! As always your help has been extremely precious!

Yours sincerely

Marco