PETSc Error: Local columns of A10 do not equal local rows of A00 with Non conforming Finite Element

General Discussion
1

Hello all,

Inspired by what I have found on the FreeFEM git repo, I am solving Stokes equation in 3D in parallel using a fieldsplit preconditioner.

  • I note my vectorial FE Pk4 = [PkV, PkV, PkV, PkM] where PkV is the FE for the velocity and PkM the FE for the pressure.

  • When PkV and PkM are classical FE let’s say P1, P2, P3, ... the code works.

  • However we have developed a nonconforming FE, called P2pnc_3d in FreeFEM, which is a non conforming finite element. Then when I use PkV=P2pnc3d and PkM=P1dc, I get the following error with PETSc:

[2]PETSC ERROR: Nonconforming object sizes
[2]PETSC ERROR: Local columns of A10 17454 do not equal local rows of A00 2291

We have tried to understand what could be this error, but we found nothing in the PETSc documention.

Could someone has already find this error ?

Besides we have another interrogation to better understand: in the code when we use Fieldsplitting, where is defined the different splits : how the solver know which block is the velocity and which block is the pressure ? The option Vh def(uh) = [1.0, 1.0, 1.0, 2.0]; and fields = uh[] let me think that there is a kind of labelling (1 for Velocity block and 2 for Pressure block), is it, or is it just a kind of initilialization ?

Thank you in advance,

Best regards,

Loïc,

Here my code:

Stokes_PAR_PETSc_3D.edp (6.7 KB)

That I run just for a small test as: mpiexec -n 4 FreeFem++-mpi Stokes_PAR_PETSc_3D.edp -split 4

2

You must define uh after (not before) you have partitioned your domain.

3

I don’t know what is up with your discretization, don’t have much time to investigate right now, but I’d suggest to:

  • remove the epsq * p * q from the varf;
  • do not pass the fields = ..., names = ... parameters in the set();
  • change the prefixes of the sub PC (from velocity_ to 0_ and pressure_ to 1_).

Then the linear system is being solved.

4

Hello @prj ,

Thanks for your reply,

It works with you solution.

I have tried to understand what this changes by studying PETSc documentation,

https://petsc.org/release/docs/manualpages/PC/PCFIELDSPLIT/
https://petsc.org/release/docs/manual/ksp/#sec-block-matrices

but I hardly found answers to my questions.

  • Why do removing epsq *p *q from the varf make the code works ?

  • Why the previous code works with the classical P2 and not for the non conforming one ? I do not see where the characteristics of the FE intervenes in the solver.

Thank you in advance for you help,

Best regards,

Loïc,

5
  • With the epsq * p * q term, your problem is a not a saddle-point problem, so the option -pc_fieldsplit_detect_saddle_point_problem does not work.
  • Your finite element space is wrong, I think. If you output the value of uh, you should only see 1 and 2 on your terminal, but in practice, you’ll get
[...]
3.3333333333e-01	3.3333333333e-01	3.3333333333e-01	3.3333333333e-01	3.3333333333e-01
	3.3333333333e-01	3.3333333333e-01	3.3333333333e-01	1.0000000000e+00	1.0000000000e+00
	1.0000000000e+00	2.0000000000e+00	2.0000000000e+00	2.0000000000e+00	
[...]

so you cannot pass this to fieldsplit because this is not a proper splitting of your fespace. There is probably something wrong with your interpolation.

6

I spoke too quickly. Because your fespace is nonconforming, there is no one-to-one correspondance between the set of ones and twos with velocity and pressure. You either need to find such a correspondance, or let PETSc detect the saddle-point structure of your discretization.

7

Thanks for your reply,

If I understand well, the labelling uh=[1.0, 1.0, 1.0, 2.0] allows to have the proper splitting of the fespace; is it ? So even if the problem is not a saddle point problem, with for example the Stokes code in github, PETSc is able to split properly the fespace.

However, when we cannot do this kind of labelling, PETSc needs to detect automatically the saddle point structure (that’s why we remove the term epspq * p * q ), is it ?
In this case, is PETSc able to determine which split is the velocity and which one is the pressure ? (and how) ?

Concerning your last remark, how can we find such correspondence ? What will be the advantage to do this ?

Best regards,

Loïc,

8

If I understand well, the labelling uh=[1.0, 1.0, 1.0, 2.0] allows to have the proper splitting of the fespace; is it ?

Yes. With a collection of conforming fespace, all fespace that have the same index are grouped in the same PETSc field.

However, when we cannot do this kind of labelling, PETSc needs to detect automatically the saddle point structure (that’s why we remove the term epspq * p * q ), is it ?

Right, here, you are “lucky” because you have only two fields, so the saddle-point structure exactly corresponds to velocity and pressure. If you want to couple Stokes with, e.g., diffusion, then you would get both velocity and diffusion in the same field, why may be sub-optimal.

Concerning your last remark, how can we find such correspondence ?

That’s for you to tell me, I don’t know how you designed the fespace.

What will be the advantage to do this ?

None if you stick to “just” Stokes.

If you are serious about preconditioning this efficiently, get in touch with me in private and we can setup a meeting.

9

Thank you for your help,

I have sent to you a mail to discuss about this,

Best regards,

Loïc,