schurPreconditioner & schurList

General Discussion
1

It seems very hard to get the the parameters schurPreconditioner and schurList correct for a PETSc block preconditioner. For a block matrix like:

image
image353Ă—151 2.93 KB

Is the following syntax correct for prameters schurPreconditioner and schurList?

func Pk=[P2,P2,P2];

fespace Sh(ThS, Pk);
fespace Vh(Th, Pk);

fespace Ph(Th, P1);

real[int] list(tndof);

list(0 : u2dof-1) = 0;
list(u2dof : u2dof+pndof - 1) = 1:pndof;
list(u2dof+pndof:u2dof+p2dof-1) = pndof+1:p2dof;
list(u2dof+p2dof:u2dof+p2dof+andof-1) = p2dof+1:p2dof+andof;
list(u2dof+p2dof+andof:tndof-1) = p2dof+andof+1 : p2dof+a2dof;

matrix[int] S(4);

varf vSchur(p, q) = int3(p * q);
S[0] = vSchur(Ph, Ph);
S[1] =S[0];

varf aSchur([ax,ay,az],[ahx,ahy,ahz]) = int2d(ThS)( vect(ah)'vect(a) );
S[2] = aSchur(Sh,Sh);
S[3] = S[2];

set(M, sparams = …
fields = usPETSc, names = names, schurPreconditioner = S, schurList = list);

Thanks a lot for your help in advance!

2

The best way to know for sure is to just try the code.

3

There is an error unfortunately:

image
image960Ă—244 13.3 KB

First, I tested an simpler case by filling the diagonal with mass matrices:

Mat M=

[[UU,  0,   UP’, 0,    UA’, 0],
 [UU,  UU,  0,   UP’,  0,   UA’],
 [UP,  0,   PP,  0,    0,   0],
 [0,   UP,  0,   PP,   0,   0],
 [UA,  0,   0,   0,    AA,  0],
 [0,   UA,  0,   0,    0,   AA]];

and set the params like:
set(M, sparams = "-ksp_monitor -ksp_type fgmres -ksp_converged_reason -ksp_initial_guess_nonzero true -ksp_rtol 1.e-5 -ksp_max_it 100 "
+ "-pc_type fieldsplit -pc_fieldsplit_type multiplicative "
+ "-fieldsplit_u_pc_type gamg "
+ "-fieldsplit_v_pc_type gamg "
+ "-fieldsplit_p_pc_type gamg "
+ "-fieldsplit_q_pc_type gamg "
+ "-fieldsplit_a_pc_type gamg "
+ “-fieldsplit_b_pc_type gamg”,
fields = usPETSc, names = names);

There is no error although the solver didn’t converge.

However, when I add

, schurPreconditioner = S, schurList = list)

the above error appears. So, it seems that my S and list are not correct?

I am not sure the format of schurPreconditioner and schurList. The only reference I can find is the following paper. But the block matrix in this paper is much simpler, it has three blocks, and S[0] is the Schur preconditioner. Regarding parameter schurList, say my Schur is a 1000x1000 matrix, should I create a schurList from 1 to 1000?

image
image656Ă—428 61.8 KB

By the way, chatGTP suggested the following, but the same error occurs.

real[int] schurList(tndof);

// u, v → not Schur
schurList(0 : u2dof-1) = 0;

// p, q → Schur block 1
schurList(u2dof : u2dof+p2dof-1) = 1;

// a, b → Schur block 2
schurList(u2dof+p2dof : tndof-1) = 2;
4

Please share a minimal code to reproduce the issue.

5

stokes_adjoint_six_fields.edp (7.1 KB)

bc1.log (2.3 MB)

fluid1.mesh (4.4 MB)

Sorry, it still need to reads a mesh and boundary data from bc1.log
I can spend more time to simply the code by generating a 3D mesh in FreeFem ….
If you remove

schurPreconditioner = S, schurList = list)

there would be no error. Thank you very much for your help with this.

6

stokes_adjoint_six_fields_simple_mesh.edp (4.7 KB)

This is a simpler version.

7

The code runs fine for me. What is your FreeFEM and PETSc versions?

8

That’s strange. FreeFEM is 4.12

9

That’s not strange at all since you are using an outdated FreeFEM version. Use the develop branch and try again, please.

10

Thank you for your help. I will update my FreeFEM.

But can I confirm that when you run

stokes_adjoint_six_fields_simple_mesh.edp (4.7 KB)

, schurPreconditioner = S, schurList = schurList);

has been included?

11

It was not in the original code, so it was not included… Now, I get an error indeed. I will look into it.

1 Like
12

OK, now I see the issue. You cannot use those parameters when the Mat has no associated domain decomposition. Since you are using a block constructor for M, it indeed has no domain decomposition attached to it. In the case, I case the simplest thing to do would be to just write a custom func and use a PCSHELL.

13

I have added an assert here: Ensure that there is a domain decomposition attached to the Mat · FreeFem/FreeFem-sources@7bb1771 · GitHub, thanks for the report.