SNES for bordered block matrix

General Discussion
23

Even with one process, I don’t get a deterministic behavior.

 finding initial point 
 initial point found. moving to continuation part 
 starting loop: 1/5. Current nu = 0.1. Predicted nu = 0.09973
  0 SNES Function norm 2.805428653365e-07 
 nu = 0.09973
  1 SNES Function norm 1.695141030744e-07 
 nu = 0.0997301
  2 SNES Function norm 6.278092520575e-13 
 nu = 0.0997308
  3 SNES Function norm 1.354552924947e-14 
Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 3
 starting loop: 2/5. Current nu = 0.0997308. Predicted nu = 0.0994608
  0 SNES Function norm 1.891297254561e-07 
 nu = 0.0994608
  1 SNES Function norm 7.815465447968e-13 
 nu = 0.0994616
  2 SNES Function norm 1.375679948496e-14 
Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 2
 starting loop: 3/5. Current nu = 0.0994616. Predicted nu = 0.0991916
  0 SNES Function norm 2.822903166950e-07 
 nu = 0.0991916
  1 SNES Function norm 1.709241490889e-07 
 nu = 0.0991917
  2 SNES Function norm 6.395819877626e-13 
 nu = 0.0991924
  3 SNES Function norm 1.349931255097e-14 
Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 3

vs.

 finding initial point 
 initial point found. moving to continuation part 
 starting loop: 1/5. Current nu = 0.1. Predicted nu = 0.09973
  0 SNES Function norm 2.805428653365e-07 
 nu = 0.09973
  1 SNES Function norm 1.695141030744e-07 
 nu = 0.0997301
  2 SNES Function norm 6.278092520575e-13 
 nu = 0.0997308
  3 SNES Function norm 1.354552924947e-14 
Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 3
 starting loop: 2/5. Current nu = 0.0997308. Predicted nu = 0.0994608
  0 SNES Function norm 2.081688271474e-03 
 nu = 0.0994608
  1 SNES Function norm 5.759118650479e-09 
 nu = 0.0994141
  2 SNES Function norm 1.409543210835e-14 
Nonlinear solve converged due to CONVERGED_FNORM_RELATIVE iterations 2
 starting loop: 3/5. Current nu = 0.0994141. Predicted nu = 0.0991441
  0 SNES Function norm 1.900550339798e-07 
 nu = 0.0991441
Nonlinear solve did not converge due to DIVERGED_LINE_SEARCH iterations 0
 starting loop: 4/5. Current nu = 0.0991441. Predicted nu = 0.0988741
  0 SNES Function norm 3.809037608138e-07

If you find where it starts to behave differently, I could try to make sense of this, but I cannot debug the full code right now, sorry.

24

Interesting. Thank you. Appreciate your help. Will follow up if I find any clues.

25

Ah, OK, I see at least one simple error. You are calling J(yq,yq), but yq is the result of a KSPSolve(), i.e., it’s in PETSc numbering. For J(.,.) you need to feed vectors in FreeFEM numbering. So, either:

  • replace KSPSolve() by ^-1 or add a call to ChangeNumbering(..., inverse = true, exchange = true), or,
  • compute the scalar product yourself with a vector in PETSc numbering (tmp = yq' * yq; mpiAllReduce(...);).
26

Thank you for pointing out that issue. This does not solve the non-determinism I’m seeing, but I think it must be related to some similar bug I have involving a mix between PETSc and FreeFEM numbering. I noticed a couple other mix-ups as well, and I think I’m getting close to a working code. I should be able to figure it out from here, but thank you so much for all your help! :slight_smile:

27

Hi there Pierre, I finally found the time to get back to this, and I got both the “simple” 2x2 block system (N+1)x(N+1) and the more complex 3x3 block system (N+N+1)x(N+N+1) to work. For once, I don’t have any problems or questions for you, and instead just wanted to report back to say thank you, and all is well! :slight_smile:

Here’s the route I took:

In the 2x2 case, it was relatively straighforward using -pc_fieldsplit_type schur as you showed me. Here we have a block structure like

Mat JlPM(JlP.n, mpirank == 0 ? 1 : 0), yqPM(yqP.n, mpirank == 0 ? 1 : 0);
Ja = [[J,JlPM],[yqPM',-1.0]];

which can be solved using

set(Ja, sparams = "-ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition self ", setup = 1);
set(J, sparams = " -fieldsplit_0_ksp_type preonly -fieldsplit_0_pc_type lu ", prefix = "fieldsplit_0_", parent = Ja);

However, the 3x3 case was a bit trickier. Eventually, I gave up on dealing with it directly as a 3x3 and did it like this:

Mat JJ = [[J, 0],[H, J]];
Mat JlPM(2*JlP.n, mpirank == 0 ? 1 : 0), qmPM(2*qmP.n, mpirank == 0 ? 1 : 0);
Mat JJa = [[JJ,JlPM],[qqmPM',0]];

with

  func real[int] blockLU(real[int]& qqm) {
      PetscScalar[int] out1, out2, in1 = qqm(0:J.n-1), in2 = qqm(J.n:2*J.n-1);
      KSPSolve(J, in1, out1);
      MatMult(H, out1, out2);
      in2 -= out2;
      KSPSolve(J, in2, out2);
      PetscScalar[int] outPETSc = [out1, out2];
      return outPETSc;
    }
set(JJa, sparams = "-ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition self ", setup = 1);
set(JJ, sparams = "-fieldsplit_0_ksp_type preonly -fieldsplit_0_pc_type shell ", precon = blockLU, prefix = "fieldsplit_0_", parent = JJa);
set(J, sparams = "-fieldsplit_ksp_type preonly -fieldsplit_pc_type lu", parent = JJ);

Thanks again!

28

Eheh, a nest of nest is something completely valid which can be useful to breakdown complexity, as seen right here. Nicely done!

1 Like