[petsc-users] Trying to apply FieldSplitPC by reading bloc matrix

Fri Dec 19 08:14:10 CST 2014

Hello Petsc Users,

I have an issue trying to use FiledSplitPC in parallel.

My goal : I want to get a linear system from petsc binary files and 
solve this in parallel with the FieldSplitPC.

The problem I want to solve is an FE approximation of the Stokes equations.

Skipping the details, my code looks like :

// Reading the four blocs UU, UP, PU and PP
for(int i=0; i < 4; ++i)
{
       string name = matrix + to_string(i) + ".petscbin";
       PetscViewer    PETSC_matreader;
       PetscViewerBinaryOpen(PETSC_COMM_WORLD, name.c_str(), 
FILE_MODE_READ, &PETSC_matreader);
       MatCreate(PETSC_COMM_WORLD,&PETSC_subA[i]);
       MatLoad(PETSC_subA[i],PETSC_matreader);
       PetscViewerDestroy(&PETSC_matreader);
}

// Reading the RHS vector and duplicating it to create the solution vector
PetscViewerBinaryOpen(PETSC_COMM_WORLD, rhs.c_str(), FILE_MODE_READ, 
&PETSC_vecreader);
VecCreate(PETSC_COMM_WORLD,&PETSC_rhs);
VecLoad(PETSC_rhs,PETSC_vecreader);
PetscViewerDestroy(&PETSC_vecreader);
VecDuplicate(PETSC_rhs,&PETSC_sol);

// Create global matrixwith MatCreateNest
MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, PETSC_subA, &PETSC_A);
MatNestGetISs(PETSC_A, PETSC_isg, NULL);

// Setting up the ksp and precond
KSPCreate(PETSC_COMM_WORLD,&PETSC_ksp);
KSPSetOperators(PETSC_ksp,PETSC_A,PETSC_A);
KSPSetFromOptions(PETSC_ksp);

KSPGetPC(PETSC_ksp, &PETSC_pc);
PCSetType(PETSC_pc, PCFIELDSPLIT);
PCFieldSplitSetIS(PETSC_pc, "0", PETSC_isg[0]);
PCFieldSplitSetIS(PETSC_pc, "1", PETSC_isg[1]);
PCSetFromOptions(PETSC_pc);

// Solving the system and writing back the solution in rhs file
KSPSolve(PETSC_ksp,PETSC_rhs,PETSC_sol);

PetscViewer    PETSC_vecwriter;
PetscViewerBinaryOpen(PETSC_COMM_WORLD, rhs.c_str(), FILE_MODE_WRITE, 
&PETSC_vecwriter);
VecView(PETSC_sol,PETSC_vecwriter);
PetscViewerDestroy(&PETSC_vecwriter);

When I run it with 1 proc, everything works fine and I get the correct 
solution :
   0 KSP preconditioned resid norm 1.271697253018e+03 true resid norm 
2.400000000000e+01 ||r(i)||/||b|| 1.000000000000e+00
   1 KSP preconditioned resid norm 5.009545069728e+01 true resid norm 
9.166803391041e-02 ||r(i)||/||b|| 3.819501412934e-03
   2 KSP preconditioned resid norm 6.460631387766e+00 true resid norm 
4.995542253831e-02 ||r(i)||/||b|| 2.081475939096e-03
   3 KSP preconditioned resid norm 1.155895209298e+00 true resid norm 
1.515734830704e-02 ||r(i)||/||b|| 6.315561794600e-04
   4 KSP preconditioned resid norm 7.407384739634e-02 true resid norm 
9.992802256200e-03 ||r(i)||/||b|| 4.163667606750e-04
   5 KSP preconditioned resid norm 1.574456882990e-02 true resid norm 
9.994876664681e-03 ||r(i)||/||b|| 4.164531943617e-04
   6 KSP preconditioned resid norm 2.383022349902e-03 true resid norm 
9.990760645581e-03 ||r(i)||/||b|| 4.162816935659e-04
   7 KSP preconditioned resid norm 6.175379834254e-04 true resid norm 
9.990821066459e-03 ||r(i)||/||b|| 4.162842111025e-04
   8 KSP preconditioned resid norm 6.867982689960e-05 true resid norm 
9.990532094790e-03 ||r(i)||/||b|| 4.162721706163e-04
   9 KSP preconditioned resid norm 1.041091257246e-05 true resid norm 
9.990558069113e-03 ||r(i)||/||b|| 4.162732528797e-04
  10 KSP preconditioned resid norm 1.447793722489e-06 true resid norm 
9.990557786778e-03 ||r(i)||/||b|| 4.162732411158e-04
  11 KSP preconditioned resid norm 2.139317335854e-07 true resid norm 
9.990557262754e-03 ||r(i)||/||b|| 4.162732192814e-04
  12 KSP preconditioned resid norm 4.383129810322e-08 true resid norm 
9.990557306920e-03 ||r(i)||/||b|| 4.162732211217e-04
  13 KSP preconditioned resid norm 3.351461304399e-09 true resid norm 
9.990557311707e-03 ||r(i)||/||b|| 4.162732213211e-04
  14 KSP preconditioned resid norm 5.169032607321e-10 true resid norm 
9.990557312817e-03 ||r(i)||/||b|| 4.162732213674e-04

[14:49:10::INFO   ] System Solved. Final tolerance reached is 
5.16903e-10 in 14 iterations.

But if I do it with 2 procs, the resolution seems fine but the solution 
is wrong :
   0 KSP preconditioned resid norm 1.247694088756e+03 true resid norm 
2.400000000000e+01 ||r(i)||/||b|| 1.000000000000e+00
   1 KSP preconditioned resid norm 4.481954484303e+01 true resid norm 
5.277507840772e-01 ||r(i)||/||b|| 2.198961600321e-02
   2 KSP preconditioned resid norm 1.110647693456e+01 true resid norm 
4.005558168981e-02 ||r(i)||/||b|| 1.668982570409e-03
   3 KSP preconditioned resid norm 1.220368027409e+00 true resid norm 
1.877650834971e-02 ||r(i)||/||b|| 7.823545145714e-04
   4 KSP preconditioned resid norm 2.834261749922e-01 true resid norm 
1.613967205264e-02 ||r(i)||/||b|| 6.724863355265e-04
   5 KSP preconditioned resid norm 4.215090288154e-02 true resid norm 
1.562561614611e-02 ||r(i)||/||b|| 6.510673394212e-04
   6 KSP preconditioned resid norm 1.209476134754e-02 true resid norm 
1.563808960492e-02 ||r(i)||/||b|| 6.515870668718e-04
   7 KSP preconditioned resid norm 2.038835108629e-03 true resid norm 
1.564163643064e-02 ||r(i)||/||b|| 6.517348512765e-04
   8 KSP preconditioned resid norm 1.928844666836e-04 true resid norm 
1.564072761376e-02 ||r(i)||/||b|| 6.516969839065e-04
   9 KSP preconditioned resid norm 3.138911950605e-05 true resid norm 
1.564047323377e-02 ||r(i)||/||b|| 6.516863847403e-04
  10 KSP preconditioned resid norm 4.950062975470e-06 true resid norm 
1.564048216528e-02 ||r(i)||/||b|| 6.516867568865e-04
  11 KSP preconditioned resid norm 7.677242244159e-07 true resid norm 
1.564049253364e-02 ||r(i)||/||b|| 6.516871889019e-04
  12 KSP preconditioned resid norm 1.870521888617e-07 true resid norm 
1.564049269566e-02 ||r(i)||/||b|| 6.516871956526e-04
  13 KSP preconditioned resid norm 3.077235724319e-08 true resid norm 
1.564049264800e-02 ||r(i)||/||b|| 6.516871936666e-04
  14 KSP preconditioned resid norm 6.584409191524e-09 true resid norm 
1.564049264183e-02 ||r(i)||/||b|| 6.516871934095e-04
  15 KSP preconditioned resid norm 1.091619359913e-09 true resid norm 
1.564049263170e-02 ||r(i)||/||b|| 6.516871929874e-04

[15:10:58::INFO   ] System Solved. Final tolerance reached is 
1.09162e-09 in 15 iterations.

Any idea of what is wrong with this ? Is it the code or the base concept ?

Thank you.
Gilles