ChangeNumbering with nested PETSc matrix returns a zero matrix

General Discussion
1

Dear FreeFEM community,

I am encountering what looks to be a bug when using ChangeNumbering with a nested block matrix corresponding to the discretization of the stokes system.

In the example attached, I create a Stokes PETSc matrix A and I supply a matrix of right hand-sides RHS2. I display the infinity norm. Then I call ChangeNumbering to convert this matrix to a new one RHSB and I observe that the infinity norm is zero.
Trying them to solve the Stokes system with multiple right hand sides then returns zero, due to this zero right hand side. What is weird is that the problem arises only if the matrix with multiple right-hand side is of dimension larger than 4.

Did I miss something?

disp("RHS2.linfty="+RHS2.linfty); //Should show 1
ChangeNumbering(A, RHS2, RHSB); 
disp("RHSB.linfty="+RHSB.linfty); //Should show 1 but shows zero
set(A, sparams = "-ksp_type preonly -pc_type lu -ksp_monitor -ksp_view -pc_mat_solver_type mumps");
KSPSolve(A,RHSB,X);  //Returns zero

The file is there:
example.zip (513.1 KB)

You should run it with

ff-mpirun -np 1 run/stokes.edp

Thanks a lot for your help.

2

Edit: I found a solution by calling ChangeNumbering on the subarrays corresponding to the block submatrices of the stokes system. I can then solve the linear system with multiple right-hand sides and call ChangeNumbering back on the sub arrays of the solution.
So I guess ChangeNumbering does not support nested block matrices yet.

matrixToIntInt(rhs2, RHS2); 
    
real[int,int] RHS2Petsc(AAstokes.n, RHS2.m);    
ChangeNumbering(AAstokes, RHS2, RHS2Petsc); 
real[int,int] RHSPPetsc(Fh1.ndof+1,RHS2.m); 
    
matrix RHSPetsc = [[RHS2Petsc],  
                   [RHSPPetsc]]; 
real[int,int] rhsPetsc(Fhv.ndof+Fh1.ndof+1,RHS2.m);     
matrixToIntInt(RHSPetsc, rhsPetsc); 

real[int,int] solV(Fhv.ndof, RHS2.m);   
real[int,int] solP(Fh1.ndof,RHS2.m); 
real[int,int] X(A.n,RHS2.m); 
//cout << "RHSB=" << RHSB << endl; 
dispVar(rhsPetsc.linfty); 
set(A, sparams = "-ksp_type preonly -pc_type lu -ksp_monitor -ksp_view -pc_mat_solver_type mumps");
KSPSolve(A,rhsPetsc,X); 
    
real[int,int] solVPetsc(AAstokes.n,RHS2.m), solPPetsc(dC.n,RHS2.m); 
solVPetsc = X(0:AAstokes.n-1,0:RHS2.m-1); 
solPPetsc = X(AAstokes.n:X.n-1,0:RHS2.m-1); 

ChangeNumbering(AAstokes, solV, solVPetsc, inverse = true, exchange = true);
ChangeNumbering(dC, solP, solPPetsc, inverse = true, exchange = true);
matrix SOLV = solV;     
matrix SOLP = solP;
3

Two comments:

  1. it is not possible to use ChangeNumbering() with a MatNest for which all diagonal blocks are not coming from a MatCreate() (I think this answers your question as to why your code what initially failing since you have 0s on the diagonal);
  2. the ChangeNumbering() has been designed to make the conversion easier for coupled problems, you should be able to do ChangeNumbering([AAstokes, dC], [solV, solP], X, inverse = true, exchange = true);, see, e.g., FreeFem-sources/examples/hpddm/stokes-block-2d-PETSc.edp at master · FreeFem/FreeFem-sources · GitHub instead of having a single call to ChangeNumbering() for each block.
4

Thanks for your reply and the explanation @prj !
I tried your solution but it seems it does not work if solV and solP are of types real[int,int].
I get this error:

131 : ChangeNumbering([AAstokes, dC], [solV, solP], X, inverse = true, exchange = true) error operator (  <7E_Array>, <7E_Array>, <P3KNMIdE> 
 List of choices 
	 (	  <l> :   <N5PETSc14DistributedCSRIN5HPDDM7SchwarzIdEEEE>, <P2KNIdE>, <P2KNIdE> )
	 (	  <l> :   <N5PETSc14DistributedCSRIN5HPDDM7SchwarzIdEEEE>, <13FEbaseArrayKnIdE>, <P3KNMIdE> )
	 (	  <l> :   <N5PETSc14DistributedCSRIN5HPDDM7SchwarzIdEEEE>, <P3KNMIdE>, <P3KNMIdE> )
	 (	  <l> :   <N5PETSc14DistributedCSRIN5HPDDM7SchwarzIdEEEE>, <P2KNIdE>, <P2KNIlE> )
	 (	  <l> :   <7E_Array>, <7E_Array>, <P2KNIdE> )
	 (	  <l> :   <P2KNIlE>, <P2KNIdE>, <P2KNIdE> )
 Error line number 131, in file run/stokes.edp, before  token )
5

Right, they must be of single dimension, I didn’t catch that you were using two-dimensional arrays. I can try to have a look to see how difficult it would be to add support for two-dimensional arrays.

6

Are you actually using this code? There are multiple errors, worth addressing if you want to use it.

7

Dear @prj, yes I am using it.
Please find attached a revised version. What errors do you see?

Thanks a lot for the help.
Regards
example.zip (66.3 KB)

8

The file run/ffimport/matrix_rhs2 is missing from your latest .zip. I see at least two issues:

  1. you are not initializing RHSPPetsc, so there may be nan in the RHS, you should add an explicit RHSPPetsc = 0
  2. solPPetsc = X(AAstokes.n:X.n-1,0:RHS2.m-1); should read solPPetsc = X(AAstokes.n:X.n-2,0:RHS2.m-1);
    If you send me an updated (and running) code, I can make sure the new ChangeNumbering() (with real[int, int] is functioning properly – the code is there but it’s not pushed yet).
9

Dear @prj, thanks a lot, indeed these were mistakes, and this resolved a bug in my application.
I was confused with the initialization because the syntax

real[int] RHSPPetsc(Fhp.ndof) = 0;

is invalid.

Here is an updated running code:
example.zip (3.0 MB) (run with ff-mpirun -np 1 run/stokes.edp)

Thanks a lot for your assistance.

10

I don’t think your code runs with more than one process, so I tested the new code in some other file, it should be good to go, you’ll need Add support for ChangeNumbering with array of Mat and KNM · FreeFem/FreeFem-sources@0c1774d · GitHub.

11

Yes indeed, I don’t need more than one process for the moment.
Thank you for the update of FreeFEM sources!