Thank you Barry for the quick reply. About the null space: I already tried what you suggest, building 2 Vec (constants) with 0 and 1 chosen by sub-domain, normalizing them and setting the null space like this MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,nconstants,constants,&nullspace); The solution is slightly different in values but it is still different in the two sub-domains. About the solver: I tried BCGS, GMRES and FGMRES. The linear system is a pressure system in a navier-stokes solver and only solving with FGMRES makes the CFD stable, with BCGS and GMRES the CFD solution diverges. Moreover, in the same case but with a single domain, CFD solution is stable using all the solvers, but FGMRES converges in much less iterations than the others. Marco Cisternino From: Barry Smith <[email protected]> Sent: mercoledì 29 settembre 2021 15:59 To: Marco Cisternino <[email protected]> Cc: [email protected] Subject: Re: [petsc-users] Disconnected domains and Poisson equation The problem actually has a two dimensional null space; constant on each domain but possibly different constants. I think you need to build the MatNullSpace by explicitly constructing two vectors, one with 0 on one domain and constant value on the other and one with 0 on the other domain and constant on the first. Separate note: why use FGMRES instead of just GMRES? If the problem is linear and the preconditioner is linear (no GMRES inside the smoother) then you can just use GMRES and it will save a little space/work and be conceptually clearer. Barry On Sep 29, 2021, at 8:46 AM, Marco Cisternino <[email protected]<mailto:[email protected]>> wrote: Good morning, I want to solve the Poisson equation on a 3D domain with 2 non-connected sub-domains. I am using FGMRES+GAMG and I have no problem if the two sub-domains see a Dirichlet boundary condition each. On the same domain I would like to solve the Poisson equation imposing periodic boundary condition in one direction and homogenous Neumann boundary conditions in the other two directions. The two sub-domains are symmetric with respect to the separation between them and the operator discretization and the right hand side are symmetric as well. It would be nice to have the same solution in both the sub-domains. Setting the null space to the constant, the solver converges to a solution having the same gradients in both sub-domains but different values. Am I doing some wrong with the null space? I’m not setting a block matrix (one block for each sub-domain), should I? I tested the null space against the matrix using MatNullSpaceTest and the answer is true. Can I do something more to have a symmetric solution as outcome of the solver? Thank you in advance for any comments and hints. Best regards, Marco Cisternino

On Sep 29, 2021, at 11:59 AM, Marco Cisternino <[email protected]> wrote:

For sake of completeness, explicitly building the null space using a vector per sub-domain make s the CFD runs using BCGS and GMRES more stable, but still slower than FGMRES.

Something is strange. Please run with -ksp_view and send the output on the solver details.

I had divergence using BCGS and GMRES setting the null space with only one constant. Thanks

Marco Cisternino

From: Marco Cisternino Sent: mercoledì 29 settembre 2021 17:54 To: Barry Smith <[email protected]> Cc: [email protected] Subject: RE: [petsc-users] Disconnected domains and Poisson equation

Thank you Barry for the quick reply. About the null space: I already tried what you suggest, building 2 Vec (constants) with 0 and 1 chosen by sub-domain, normalizing them and setting the null space like this MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,nconstants,constants,&nullspace); The solution is slightly different in values but it is still different in the two sub-domains. About the solver: I tried BCGS, GMRES and FGMRES. The linear system is a pressure system in a navier-stokes solver and only solving with FGMRES makes the CFD stable, with BCGS and GMRES the CFD solution diverges. Moreover, in the same case but with a single domain, CFD solution is stable using all the solvers, but FGMRES converges in much less iterations than the others.

Marco Cisternino

From: Barry Smith <[email protected] <mailto:[email protected]>> Sent: mercoledì 29 settembre 2021 15:59 To: Marco Cisternino <[email protected] <mailto:[email protected]>> Cc: [email protected] <mailto:[email protected]> Subject: Re: [petsc-users] Disconnected domains and Poisson equation

The problem actually has a two dimensional null space; constant on each domain but possibly different constants. I think you need to build the MatNullSpace by explicitly constructing two vectors, one with 0 on one domain and constant value on the other and one with 0 on the other domain and constant on the first.

Separate note: why use FGMRES instead of just GMRES? If the problem is linear and the preconditioner is linear (no GMRES inside the smoother) then you can just use GMRES and it will save a little space/work and be conceptually clearer.

Barry

On Sep 29, 2021, at 8:46 AM, Marco Cisternino <[email protected] <mailto:[email protected]>> wrote:

Good morning, I want to solve the Poisson equation on a 3D domain with 2 non-connected sub-domains. I am using FGMRES+GAMG and I have no problem if the two sub-domains see a Dirichlet boundary condition each. On the same domain I would like to solve the Poisson equation imposing periodic boundary condition in one direction and homogenous Neumann boundary conditions in the other two directions. The two sub-domains are symmetric with respect to the separation between them and the operator discretization and the right hand side are symmetric as well. It would be nice to have the same solution in both the sub-domains. Setting the null space to the constant, the solver converges to a solution having the same gradients in both sub-domains but different values. Am I doing some wrong with the null space? I’m not setting a block matrix (one block for each sub-domain), should I? I tested the null space against the matrix using MatNullSpaceTest and the answer is true. Can I do something more to have a symmetric solution as outcome of the solver? Thank you in advance for any comments and hints.

Best regards,

Marco Cisternino

* Do we understand: type: chebyshev eigenvalue estimates used: *min = 0., max = 0.* eigenvalues estimate via gmres *min 0., max 0.* * Is this Poisson solver unsymmetric? * Does this problem start off converging and then evolve and then start stagnating? The eigen estimates may need to be recalculated. Also Chebyshev is problematic for unsymmetric matrices. Hypre does better with unsymmetric matrices. Mark On Thu, Sep 30, 2021 at 4:16 AM Marco Cisternino < [email protected]> wrote:

Hello Barry.

This is the output of ksp_view using fgmres and gamg. It has to be said that the solution of the linear system should be a zero values field. As you can see both unpreconditioned residual and r/b converge at this iteration of the CFD solver. During the time integration of the CFD, I can observe pressure linear solver residuals behaving in a different way: unpreconditioned residual stil converges but r/b stalls. After the output of ksp_view I add the output of ksp_monitor_true_residual for one of these iteration where r/b stalls. Thanks,

KSP Object: 1 MPI processes

type: fgmres

restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement

happy breakdown tolerance 1e-30

maximum iterations=100, nonzero initial guess

tolerances: relative=1e-05, absolute=1e-50, divergence=10000.

right preconditioning

using UNPRECONDITIONED norm type for convergence test

PC Object: 1 MPI processes

type: gamg

type is MULTIPLICATIVE, levels=4 cycles=v

Cycles per PCApply=1

Using externally compute Galerkin coarse grid matrices

GAMG specific options

Threshold for dropping small values in graph on each level = 0.02 0.02

Threshold scaling factor for each level not specified = 1.

AGG specific options

Symmetric graph true

Number of levels to square graph 1

Number smoothing steps 0

Coarse grid solver -- level -------------------------------

KSP Object: (mg_coarse_) 1 MPI processes

type: preonly

maximum iterations=10000, initial guess is zero

tolerances: relative=1e-05, absolute=1e-50, divergence=10000.

left preconditioning

using NONE norm type for convergence test

PC Object: (mg_coarse_) 1 MPI processes

type: bjacobi

number of blocks = 1

Local solve is same for all blocks, in the following KSP and PC objects:

KSP Object: (mg_coarse_sub_) 1 MPI processes

type: preonly

maximum iterations=1, initial guess is zero

tolerances: relative=1e-05, absolute=1e-50, divergence=10000.

left preconditioning

using DEFAULT norm type for convergence test

PC Object: (mg_coarse_sub_) 1 MPI processes

type: lu

PC has not been set up so information may be incomplete

out-of-place factorization

tolerance for zero pivot 2.22045e-14

using diagonal shift on blocks to prevent zero pivot [INBLOCKS]

matrix ordering: nd

linear system matrix = precond matrix:

Mat Object: 1 MPI processes

type: seqaij

rows=18, cols=18

total: nonzeros=104, allocated nonzeros=104

total number of mallocs used during MatSetValues calls =0

not using I-node routines

linear system matrix = precond matrix:

Mat Object: 1 MPI processes

type: seqaij

rows=18, cols=18

total: nonzeros=104, allocated nonzeros=104

total number of mallocs used during MatSetValues calls =0

not using I-node routines

Down solver (pre-smoother) on level 1 -------------------------------

KSP Object: (mg_levels_1_) 1 MPI processes

type: chebyshev

eigenvalue estimates used: min = 0., max = 0.

eigenvalues estimate via gmres min 0., max 0.

eigenvalues estimated using gmres with translations [0. 0.1; 0. 1.1]

KSP Object: (mg_levels_1_esteig_) 1 MPI processes

type: gmres

restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement

happy breakdown tolerance 1e-30

maximum iterations=10, initial guess is zero

tolerances: relative=1e-12, absolute=1e-50, divergence=10000.

left preconditioning

using DEFAULT norm type for convergence test

estimating eigenvalues using noisy right hand side

maximum iterations=2, nonzero initial guess

tolerances: relative=1e-05, absolute=1e-50, divergence=10000.

left preconditioning

using NONE norm type for convergence test

PC Object: (mg_levels_1_) 1 MPI processes

type: sor

type = local_symmetric, iterations = 1, local iterations = 1, omega = 1.

linear system matrix = precond matrix:

Mat Object: 1 MPI processes

type: seqaij

rows=67, cols=67

total: nonzeros=675, allocated nonzeros=675

total number of mallocs used during MatSetValues calls =0

not using I-node routines

Up solver (post-smoother) same as down solver (pre-smoother)

Down solver (pre-smoother) on level 2 -------------------------------

KSP Object: (mg_levels_2_) 1 MPI processes

type: chebyshev

eigenvalue estimates used: min = 0., max = 0.

eigenvalues estimate via gmres min 0., max 0.

eigenvalues estimated using gmres with translations [0. 0.1; 0. 1.1]

KSP Object: (mg_levels_2_esteig_) 1 MPI processes

type: gmres

restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement

happy breakdown tolerance 1e-30

maximum iterations=10, initial guess is zero

tolerances: relative=1e-12, absolute=1e-50, divergence=10000.

left preconditioning

using DEFAULT norm type for convergence test

estimating eigenvalues using noisy right hand side

maximum iterations=2, nonzero initial guess

tolerances: relative=1e-05, absolute=1e-50, divergence=10000.

left preconditioning

using NONE norm type for convergence test

PC Object: (mg_levels_2_) 1 MPI processes

type: sor

type = local_symmetric, iterations = 1, local iterations = 1, omega = 1.

linear system matrix = precond matrix:

Mat Object: 1 MPI processes

type: seqaij

rows=348, cols=348

total: nonzeros=3928, allocated nonzeros=3928

total number of mallocs used during MatSetValues calls =0

not using I-node routines

Up solver (post-smoother) same as down solver (pre-smoother)

Down solver (pre-smoother) on level 3 -------------------------------

KSP Object: (mg_levels_3_) 1 MPI processes

type: chebyshev

eigenvalue estimates used: min = 0., max = 0.

eigenvalues estimate via gmres min 0., max 0.

eigenvalues estimated using gmres with translations [0. 0.1; 0. 1.1]

KSP Object: (mg_levels_3_esteig_) 1 MPI processes

type: gmres

restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement

happy breakdown tolerance 1e-30

maximum iterations=10, initial guess is zero

tolerances: relative=1e-12, absolute=1e-50, divergence=10000.

left preconditioning

using DEFAULT norm type for convergence test

estimating eigenvalues using noisy right hand side

maximum iterations=2, nonzero initial guess

tolerances: relative=1e-05, absolute=1e-50, divergence=10000.

left preconditioning

using NONE norm type for convergence test

PC Object: (mg_levels_3_) 1 MPI processes

type: sor

type = local_symmetric, iterations = 1, local iterations = 1, omega = 1.

linear system matrix = precond matrix:

Mat Object: 1 MPI processes

type: seqaij

rows=3584, cols=3584

total: nonzeros=23616, allocated nonzeros=23616

total number of mallocs used during MatSetValues calls =0

has attached null space

not using I-node routines

Up solver (post-smoother) same as down solver (pre-smoother)

linear system matrix = precond matrix:

Mat Object: 1 MPI processes

type: seqaij

rows=3584, cols=3584

total: nonzeros=23616, allocated nonzeros=23616

total number of mallocs used during MatSetValues calls =0

has attached null space

not using I-node routines

Pressure system has reached convergence in 0 iterations with reason 3.

0 KSP unpreconditioned resid norm 4.798763170703e-16 true resid norm 4.798763170703e-16 ||r(i)||/||b|| 1.000000000000e+00

0 KSP Residual norm 4.798763170703e-16 % max 1.000000000000e+00 min 1.000000000000e+00 max/min 1.000000000000e+00

1 KSP unpreconditioned resid norm 1.648749109132e-17 true resid norm 1.648749109132e-17 ||r(i)||/||b|| 3.435779284125e-02

1 KSP Residual norm 1.648749109132e-17 % max 9.561792537103e-01 min 9.561792537103e-01 max/min 1.000000000000e+00

2 KSP unpreconditioned resid norm 4.737880600040e-19 true resid norm 4.737880600040e-19 ||r(i)||/||b|| 9.873128619820e-04

2 KSP Residual norm 4.737880600040e-19 % max 9.828636644296e-01 min 9.293131521763e-01 max/min 1.057623753767e+00

3 KSP unpreconditioned resid norm 2.542212716830e-20 true resid norm 2.542212716830e-20 ||r(i)||/||b|| 5.297641551371e-05

3 KSP Residual norm 2.542212716830e-20 % max 9.933572357920e-01 min 9.158303248850e-01 max/min 1.084652046127e+00

4 KSP unpreconditioned resid norm 6.614510286263e-21 true resid norm 6.614510286269e-21 ||r(i)||/||b|| 1.378378146822e-05

4 KSP Residual norm 6.614510286263e-21 % max 9.950912550705e-01 min 6.296575800237e-01 max/min 1.580368896747e+00

5 KSP unpreconditioned resid norm 1.981505525281e-22 true resid norm 1.981505525272e-22 ||r(i)||/||b|| 4.129200493513e-07

5 KSP Residual norm 1.981505525281e-22 % max 9.984097962703e-01 min 5.316259535293e-01 max/min 1.878030577029e+00

Linear solve converged due to CONVERGED_RTOL iterations 5

Ksp_monitor_true_residual output for stalling r/b CFD iteration 0 KSP unpreconditioned resid norm 9.010260489109e-14 true resid norm 9.010260489109e-14 ||r(i)||/||b|| 2.021559024868e+00

0 KSP Residual norm 9.010260489109e-14 % max 1.000000000000e+00 min 1.000000000000e+00 max/min 1.000000000000e+00

1 KSP unpreconditioned resid norm 4.918108339808e-15 true resid norm 4.918171792537e-15 ||r(i)||/||b|| 1.103450292594e-01

1 KSP Residual norm 4.918108339808e-15 % max 9.566256813737e-01 min 9.566256813737e-01 max/min 1.000000000000e+00

2 KSP unpreconditioned resid norm 1.443599554690e-15 true resid norm 1.444867143493e-15 ||r(i)||/||b|| 3.241731154382e-02

2 KSP Residual norm 1.443599554690e-15 % max 9.614019380614e-01 min 7.360950481750e-01 max/min 1.306083963538e+00

3 KSP unpreconditioned resid norm 6.623206616803e-16 true resid norm 6.654132553541e-16 ||r(i)||/||b|| 1.492933720678e-02

3 KSP Residual norm 6.623206616803e-16 % max 9.764112945239e-01 min 4.911485418014e-01 max/min 1.988016274960e+00

4 KSP unpreconditioned resid norm 6.551896936698e-16 true resid norm 6.646157296305e-16 ||r(i)||/||b|| 1.491144376933e-02

4 KSP Residual norm 6.551896936698e-16 % max 9.883425885532e-01 min 1.461270778833e-01 max/min 6.763582786091e+00

5 KSP unpreconditioned resid norm 6.222297644887e-16 true resid norm 1.720560536914e-15 ||r(i)||/||b|| 3.860282047823e-02

5 KSP Residual norm 6.222297644887e-16 % max 1.000409371755e+00 min 4.989767363560e-03 max/min 2.004921870829e+02

6 KSP unpreconditioned resid norm 6.496945794974e-17 true resid norm 2.031914800253e-14 ||r(i)||/||b|| 4.558842341106e-01

6 KSP Residual norm 6.496945794974e-17 % max 1.004914985753e+00 min 1.459258738706e-03 max/min 6.886475709192e+02

7 KSP unpreconditioned resid norm 1.965237342540e-17 true resid norm 1.684522207337e-14 ||r(i)||/||b|| 3.779425772373e-01

7 KSP Residual norm 1.965237342540e-17 % max 1.005737762541e+00 min 1.452603803766e-03 max/min 6.923689446035e+02

8 KSP unpreconditioned resid norm 1.627718951285e-17 true resid norm 1.958642967520e-14 ||r(i)||/||b|| 4.394448276241e-01

8 KSP Residual norm 1.627718951285e-17 % max 1.006364278765e+00 min 1.452081813014e-03 max/min 6.930492963590e+02

9 KSP unpreconditioned resid norm 1.616577677764e-17 true resid norm 2.019110946644e-14 ||r(i)||/||b|| 4.530115373837e-01

9 KSP Residual norm 1.616577677764e-17 % max 1.006648747131e+00 min 1.452031376577e-03 max/min 6.932692801059e+02

10 KSP unpreconditioned resid norm 1.285788988203e-17 true resid norm 2.065082694477e-14 ||r(i)||/||b|| 4.633258453698e-01

10 KSP Residual norm 1.285788988203e-17 % max 1.007469033514e+00 min 1.433291867068e-03 max/min 7.029057072477e+02

11 KSP unpreconditioned resid norm 5.490854431580e-19 true resid norm 1.798071628891e-14 ||r(i)||/||b|| 4.034187394623e-01

11 KSP Residual norm 5.490854431580e-19 % max 1.008058905554e+00 min 1.369401685301e-03 max/min 7.361309076612e+02

12 KSP unpreconditioned resid norm 1.371754802104e-20 true resid norm 1.965688920064e-14 ||r(i)||/||b|| 4.410256708163e-01

12 KSP Residual norm 1.371754802104e-20 % max 1.008409402214e+00 min 1.369243011779e-03 max/min 7.364721919624e+02

Linear solve converged due to CONVERGED_RTOL iterations 12

Marco Cisternino

*From:* Barry Smith <[email protected]> *Sent:* mercoledì 29 settembre 2021 18:34 *To:* Marco Cisternino <[email protected]> *Cc:* [email protected] *Subject:* Re: [petsc-users] Disconnected domains and Poisson equation

On Sep 29, 2021, at 11:59 AM, Marco Cisternino < [email protected]> wrote:

For sake of completeness, explicitly building the null space using a vector per sub-domain make s the CFD runs using BCGS and GMRES more stable, but still slower than FGMRES.

Something is strange. Please run with -ksp_view and send the output on the solver details.

I had divergence using BCGS and GMRES setting the null space with only one constant.

Thanks

Marco Cisternino

*From:* Marco Cisternino *Sent:* mercoledì 29 settembre 2021 17:54 *To:* Barry Smith <[email protected]> *Cc:* [email protected] *Subject:* RE: [petsc-users] Disconnected domains and Poisson equation

Thank you Barry for the quick reply.

About the null space: I already tried what you suggest, building 2 Vec (constants) with 0 and 1 chosen by sub-domain, normalizing them and setting the null space like this

MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,nconstants,constants,&nullspace);

The solution is slightly different in values but it is still different in the two sub-domains.

About the solver: I tried BCGS, GMRES and FGMRES. The linear system is a pressure system in a navier-stokes solver and only solving with FGMRES makes the CFD stable, with BCGS and GMRES the CFD solution diverges. Moreover, in the same case but with a single domain, CFD solution is stable using all the solvers, but FGMRES converges in much less iterations than the others.

Marco Cisternino

*From:* Barry Smith <[email protected]> *Sent:* mercoledì 29 settembre 2021 15:59 *To:* Marco Cisternino <[email protected]> *Cc:* [email protected] *Subject:* Re: [petsc-users] Disconnected domains and Poisson equation

The problem actually has a two dimensional null space; constant on each domain but possibly different constants. I think you need to build the MatNullSpace by explicitly constructing two vectors, one with 0 on one domain and constant value on the other and one with 0 on the other domain and constant on the first.

Separate note: why use FGMRES instead of just GMRES? If the problem is linear and the preconditioner is linear (no GMRES inside the smoother) then you can just use GMRES and it will save a little space/work and be conceptually clearer.

Barry

On Sep 29, 2021, at 8:46 AM, Marco Cisternino <[email protected]> wrote:

Good morning,

I want to solve the Poisson equation on a 3D domain with 2 non-connected sub-domains.

I am using FGMRES+GAMG and I have no problem if the two sub-domains see a Dirichlet boundary condition each.

On the same domain I would like to solve the Poisson equation imposing periodic boundary condition in one direction and homogenous Neumann boundary conditions in the other two directions. The two sub-domains are symmetric with respect to the separation between them and the operator discretization and the right hand side are symmetric as well. It would be nice to have the same solution in both the sub-domains.

Setting the null space to the constant, the solver converges to a solution having the same gradients in both sub-domains but different values.

Am I doing some wrong with the null space? I’m not setting a block matrix (one block for each sub-domain), should I?

I tested the null space against the matrix using MatNullSpaceTest and the answer is true. Can I do something more to have a symmetric solution as outcome of the solver?

Thank you in advance for any comments and hints.

Best regards,

Marco Cisternino

Thank you Barry. I added a custom atoll = 1.0e-12 and this makes the CFD stable with all the linear solver types. CFD solution is good and pressure is a good “zero” field at every CFD iteration. I did the same test using ASM+ILU+FGMRES(BCGS and GMRES) and the behaviour is the same. During some CFD iteration the residual of linear system starts slightly higher than atol and the linear solver makes some iteration (2/3 iterations) before it stops because of atol. The pressure is still different in the 2 sub-domains (order 1.0e-14 because of those few linear solver iterations), therefore no symmetry of the solution In the 2 sub-domains. I think it is a matter of round-off, do you agree on this? Or do I need to take care of this difference as a symptom of something wrong? Thank you for your support. Marco Cisternino From: Barry Smith <[email protected]> Sent: giovedì 30 settembre 2021 16:39 To: Marco Cisternino <[email protected]> Cc: [email protected] Subject: Re: [petsc-users] Disconnected domains and Poisson equation It looks like the initial solution (guess) is to round-off the solution to the linear system 9.010260489109e-14 0 KSP unpreconditioned resid norm 9.010260489109e-14 true resid norm 9.010260489109e-14 ||r(i)||/||b|| 2.021559024868e+00 0 KSP Residual norm 9.010260489109e-14 % max 1.000000000000e+00 min 1.000000000000e+00 max/min 1.000000000000e+00 1 KSP unpreconditioned resid norm 4.918108339808e-15 true resid norm 4.918171792537e-15 ||r(i)||/||b|| 1.103450292594e-01 1 KSP Residual norm 4.918108339808e-15 % max 9.566256813737e-01 min 9.566256813737e-01 max/min 1.000000000000e+00 2 KSP unpreconditioned resid norm 1.443599554690e-15 true resid norm 1.444867143493e-15 ||r(i)||/||b|| 3.241731154382e-02 2 KSP Residual norm 1.443599554690e-15 % max 9.614019380614e-01 min 7.360950481750e-01 max/min 1.306083963538e+00 Thus the Krylov solver will not be able to improve the solution, it then gets stuck trying to improve the solution but cannot because of round off. In other words the algorithm has converged (even at the initial solution (guess) and should stop immediately. You can use -ksp_atol 1.e-12 to get it to stop immediately without iterating if the initial residual is less than 1e-12. Barry On Sep 30, 2021, at 4:16 AM, Marco Cisternino <[email protected]<mailto:[email protected]>> wrote: Hello Barry. This is the output of ksp_view using fgmres and gamg. It has to be said that the solution of the linear system should be a zero values field. As you can see both unpreconditioned residual and r/b converge at this iteration of the CFD solver. During the time integration of the CFD, I can observe pressure linear solver residuals behaving in a different way: unpreconditioned residual stil converges but r/b stalls. After the output of ksp_view I add the output of ksp_monitor_true_residual for one of these iteration where r/b stalls. Thanks, KSP Object: 1 MPI processes type: fgmres restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement happy breakdown tolerance 1e-30 maximum iterations=100, nonzero initial guess tolerances: relative=1e-05, absolute=1e-50, divergence=10000. right preconditioning using UNPRECONDITIONED norm type for convergence test PC Object: 1 MPI processes type: gamg type is MULTIPLICATIVE, levels=4 cycles=v Cycles per PCApply=1 Using externally compute Galerkin coarse grid matrices GAMG specific options Threshold for dropping small values in graph on each level = 0.02 0.02 Threshold scaling factor for each level not specified = 1. AGG specific options Symmetric graph true Number of levels to square graph 1 Number smoothing steps 0 Coarse grid solver -- level ------------------------------- KSP Object: (mg_coarse_) 1 MPI processes type: preonly maximum iterations=10000, initial guess is zero tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using NONE norm type for convergence test PC Object: (mg_coarse_) 1 MPI processes type: bjacobi number of blocks = 1 Local solve is same for all blocks, in the following KSP and PC objects: KSP Object: (mg_coarse_sub_) 1 MPI processes type: preonly maximum iterations=1, initial guess is zero tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using DEFAULT norm type for convergence test PC Object: (mg_coarse_sub_) 1 MPI processes type: lu PC has not been set up so information may be incomplete out-of-place factorization tolerance for zero pivot 2.22045e-14 using diagonal shift on blocks to prevent zero pivot [INBLOCKS] matrix ordering: nd linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=18, cols=18 total: nonzeros=104, allocated nonzeros=104 total number of mallocs used during MatSetValues calls =0 not using I-node routines linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=18, cols=18 total: nonzeros=104, allocated nonzeros=104 total number of mallocs used during MatSetValues calls =0 not using I-node routines Down solver (pre-smoother) on level 1 ------------------------------- KSP Object: (mg_levels_1_) 1 MPI processes type: chebyshev eigenvalue estimates used: min = 0., max = 0. eigenvalues estimate via gmres min 0., max 0. eigenvalues estimated using gmres with translations [0. 0.1; 0. 1.1] KSP Object: (mg_levels_1_esteig_) 1 MPI processes type: gmres restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement happy breakdown tolerance 1e-30 maximum iterations=10, initial guess is zero tolerances: relative=1e-12, absolute=1e-50, divergence=10000. left preconditioning using DEFAULT norm type for convergence test estimating eigenvalues using noisy right hand side maximum iterations=2, nonzero initial guess tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using NONE norm type for convergence test PC Object: (mg_levels_1_) 1 MPI processes type: sor type = local_symmetric, iterations = 1, local iterations = 1, omega = 1. linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=67, cols=67 total: nonzeros=675, allocated nonzeros=675 total number of mallocs used during MatSetValues calls =0 not using I-node routines Up solver (post-smoother) same as down solver (pre-smoother) Down solver (pre-smoother) on level 2 ------------------------------- KSP Object: (mg_levels_2_) 1 MPI processes type: chebyshev eigenvalue estimates used: min = 0., max = 0. eigenvalues estimate via gmres min 0., max 0. eigenvalues estimated using gmres with translations [0. 0.1; 0. 1.1] KSP Object: (mg_levels_2_esteig_) 1 MPI processes type: gmres restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement happy breakdown tolerance 1e-30 maximum iterations=10, initial guess is zero tolerances: relative=1e-12, absolute=1e-50, divergence=10000. left preconditioning using DEFAULT norm type for convergence test estimating eigenvalues using noisy right hand side maximum iterations=2, nonzero initial guess tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using NONE norm type for convergence test PC Object: (mg_levels_2_) 1 MPI processes type: sor type = local_symmetric, iterations = 1, local iterations = 1, omega = 1. linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=348, cols=348 total: nonzeros=3928, allocated nonzeros=3928 total number of mallocs used during MatSetValues calls =0 not using I-node routines Up solver (post-smoother) same as down solver (pre-smoother) Down solver (pre-smoother) on level 3 ------------------------------- KSP Object: (mg_levels_3_) 1 MPI processes type: chebyshev eigenvalue estimates used: min = 0., max = 0. eigenvalues estimate via gmres min 0., max 0. eigenvalues estimated using gmres with translations [0. 0.1; 0. 1.1] KSP Object: (mg_levels_3_esteig_) 1 MPI processes type: gmres restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement happy breakdown tolerance 1e-30 maximum iterations=10, initial guess is zero tolerances: relative=1e-12, absolute=1e-50, divergence=10000. left preconditioning using DEFAULT norm type for convergence test estimating eigenvalues using noisy right hand side maximum iterations=2, nonzero initial guess tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using NONE norm type for convergence test PC Object: (mg_levels_3_) 1 MPI processes type: sor type = local_symmetric, iterations = 1, local iterations = 1, omega = 1. linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=3584, cols=3584 total: nonzeros=23616, allocated nonzeros=23616 total number of mallocs used during MatSetValues calls =0 has attached null space not using I-node routines Up solver (post-smoother) same as down solver (pre-smoother) linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=3584, cols=3584 total: nonzeros=23616, allocated nonzeros=23616 total number of mallocs used during MatSetValues calls =0 has attached null space not using I-node routines Pressure system has reached convergence in 0 iterations with reason 3. 0 KSP unpreconditioned resid norm 4.798763170703e-16 true resid norm 4.798763170703e-16 ||r(i)||/||b|| 1.000000000000e+00 0 KSP Residual norm 4.798763170703e-16 % max 1.000000000000e+00 min 1.000000000000e+00 max/min 1.000000000000e+00 1 KSP unpreconditioned resid norm 1.648749109132e-17 true resid norm 1.648749109132e-17 ||r(i)||/||b|| 3.435779284125e-02 1 KSP Residual norm 1.648749109132e-17 % max 9.561792537103e-01 min 9.561792537103e-01 max/min 1.000000000000e+00 2 KSP unpreconditioned resid norm 4.737880600040e-19 true resid norm 4.737880600040e-19 ||r(i)||/||b|| 9.873128619820e-04 2 KSP Residual norm 4.737880600040e-19 % max 9.828636644296e-01 min 9.293131521763e-01 max/min 1.057623753767e+00 3 KSP unpreconditioned resid norm 2.542212716830e-20 true resid norm 2.542212716830e-20 ||r(i)||/||b|| 5.297641551371e-05 3 KSP Residual norm 2.542212716830e-20 % max 9.933572357920e-01 min 9.158303248850e-01 max/min 1.084652046127e+00 4 KSP unpreconditioned resid norm 6.614510286263e-21 true resid norm 6.614510286269e-21 ||r(i)||/||b|| 1.378378146822e-05 4 KSP Residual norm 6.614510286263e-21 % max 9.950912550705e-01 min 6.296575800237e-01 max/min 1.580368896747e+00 5 KSP unpreconditioned resid norm 1.981505525281e-22 true resid norm 1.981505525272e-22 ||r(i)||/||b|| 4.129200493513e-07 5 KSP Residual norm 1.981505525281e-22 % max 9.984097962703e-01 min 5.316259535293e-01 max/min 1.878030577029e+00 Linear solve converged due to CONVERGED_RTOL iterations 5 Ksp_monitor_true_residual output for stalling r/b CFD iteration 0 KSP unpreconditioned resid norm 9.010260489109e-14 true resid norm 9.010260489109e-14 ||r(i)||/||b|| 2.021559024868e+00 0 KSP Residual norm 9.010260489109e-14 % max 1.000000000000e+00 min 1.000000000000e+00 max/min 1.000000000000e+00 1 KSP unpreconditioned resid norm 4.918108339808e-15 true resid norm 4.918171792537e-15 ||r(i)||/||b|| 1.103450292594e-01 1 KSP Residual norm 4.918108339808e-15 % max 9.566256813737e-01 min 9.566256813737e-01 max/min 1.000000000000e+00 2 KSP unpreconditioned resid norm 1.443599554690e-15 true resid norm 1.444867143493e-15 ||r(i)||/||b|| 3.241731154382e-02 2 KSP Residual norm 1.443599554690e-15 % max 9.614019380614e-01 min 7.360950481750e-01 max/min 1.306083963538e+00 3 KSP unpreconditioned resid norm 6.623206616803e-16 true resid norm 6.654132553541e-16 ||r(i)||/||b|| 1.492933720678e-02 3 KSP Residual norm 6.623206616803e-16 % max 9.764112945239e-01 min 4.911485418014e-01 max/min 1.988016274960e+00 4 KSP unpreconditioned resid norm 6.551896936698e-16 true resid norm 6.646157296305e-16 ||r(i)||/||b|| 1.491144376933e-02 4 KSP Residual norm 6.551896936698e-16 % max 9.883425885532e-01 min 1.461270778833e-01 max/min 6.763582786091e+00 5 KSP unpreconditioned resid norm 6.222297644887e-16 true resid norm 1.720560536914e-15 ||r(i)||/||b|| 3.860282047823e-02 5 KSP Residual norm 6.222297644887e-16 % max 1.000409371755e+00 min 4.989767363560e-03 max/min 2.004921870829e+02 6 KSP unpreconditioned resid norm 6.496945794974e-17 true resid norm 2.031914800253e-14 ||r(i)||/||b|| 4.558842341106e-01 6 KSP Residual norm 6.496945794974e-17 % max 1.004914985753e+00 min 1.459258738706e-03 max/min 6.886475709192e+02 7 KSP unpreconditioned resid norm 1.965237342540e-17 true resid norm 1.684522207337e-14 ||r(i)||/||b|| 3.779425772373e-01 7 KSP Residual norm 1.965237342540e-17 % max 1.005737762541e+00 min 1.452603803766e-03 max/min 6.923689446035e+02 8 KSP unpreconditioned resid norm 1.627718951285e-17 true resid norm 1.958642967520e-14 ||r(i)||/||b|| 4.394448276241e-01 8 KSP Residual norm 1.627718951285e-17 % max 1.006364278765e+00 min 1.452081813014e-03 max/min 6.930492963590e+02 9 KSP unpreconditioned resid norm 1.616577677764e-17 true resid norm 2.019110946644e-14 ||r(i)||/||b|| 4.530115373837e-01 9 KSP Residual norm 1.616577677764e-17 % max 1.006648747131e+00 min 1.452031376577e-03 max/min 6.932692801059e+02 10 KSP unpreconditioned resid norm 1.285788988203e-17 true resid norm 2.065082694477e-14 ||r(i)||/||b|| 4.633258453698e-01 10 KSP Residual norm 1.285788988203e-17 % max 1.007469033514e+00 min 1.433291867068e-03 max/min 7.029057072477e+02 11 KSP unpreconditioned resid norm 5.490854431580e-19 true resid norm 1.798071628891e-14 ||r(i)||/||b|| 4.034187394623e-01 11 KSP Residual norm 5.490854431580e-19 % max 1.008058905554e+00 min 1.369401685301e-03 max/min 7.361309076612e+02 12 KSP unpreconditioned resid norm 1.371754802104e-20 true resid norm 1.965688920064e-14 ||r(i)||/||b|| 4.410256708163e-01 12 KSP Residual norm 1.371754802104e-20 % max 1.008409402214e+00 min 1.369243011779e-03 max/min 7.364721919624e+02 Linear solve converged due to CONVERGED_RTOL iterations 12 Marco Cisternino From: Barry Smith <[email protected]<mailto:[email protected]>> Sent: mercoledì 29 settembre 2021 18:34 To: Marco Cisternino <[email protected]<mailto:[email protected]>> Cc: [email protected]<mailto:[email protected]> Subject: Re: [petsc-users] Disconnected domains and Poisson equation On Sep 29, 2021, at 11:59 AM, Marco Cisternino <[email protected]<mailto:[email protected]>> wrote: For sake of completeness, explicitly building the null space using a vector per sub-domain make s the CFD runs using BCGS and GMRES more stable, but still slower than FGMRES. Something is strange. Please run with -ksp_view and send the output on the solver details. I had divergence using BCGS and GMRES setting the null space with only one constant. Thanks Marco Cisternino From: Marco Cisternino Sent: mercoledì 29 settembre 2021 17:54 To: Barry Smith <[email protected]<mailto:[email protected]>> Cc: [email protected]<mailto:[email protected]> Subject: RE: [petsc-users] Disconnected domains and Poisson equation Thank you Barry for the quick reply. About the null space: I already tried what you suggest, building 2 Vec (constants) with 0 and 1 chosen by sub-domain, normalizing them and setting the null space like this MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,nconstants,constants,&nullspace); The solution is slightly different in values but it is still different in the two sub-domains. About the solver: I tried BCGS, GMRES and FGMRES. The linear system is a pressure system in a navier-stokes solver and only solving with FGMRES makes the CFD stable, with BCGS and GMRES the CFD solution diverges. Moreover, in the same case but with a single domain, CFD solution is stable using all the solvers, but FGMRES converges in much less iterations than the others. Marco Cisternino From: Barry Smith <[email protected]<mailto:[email protected]>> Sent: mercoledì 29 settembre 2021 15:59 To: Marco Cisternino <[email protected]<mailto:[email protected]>> Cc: [email protected]<mailto:[email protected]> Subject: Re: [petsc-users] Disconnected domains and Poisson equation The problem actually has a two dimensional null space; constant on each domain but possibly different constants. I think you need to build the MatNullSpace by explicitly constructing two vectors, one with 0 on one domain and constant value on the other and one with 0 on the other domain and constant on the first. Separate note: why use FGMRES instead of just GMRES? If the problem is linear and the preconditioner is linear (no GMRES inside the smoother) then you can just use GMRES and it will save a little space/work and be conceptually clearer. Barry On Sep 29, 2021, at 8:46 AM, Marco Cisternino <[email protected]<mailto:[email protected]>> wrote: Good morning, I want to solve the Poisson equation on a 3D domain with 2 non-connected sub-domains. I am using FGMRES+GAMG and I have no problem if the two sub-domains see a Dirichlet boundary condition each. On the same domain I would like to solve the Poisson equation imposing periodic boundary condition in one direction and homogenous Neumann boundary conditions in the other two directions. The two sub-domains are symmetric with respect to the separation between them and the operator discretization and the right hand side are symmetric as well. It would be nice to have the same solution in both the sub-domains. Setting the null space to the constant, the solver converges to a solution having the same gradients in both sub-domains but different values. Am I doing some wrong with the null space? I’m not setting a block matrix (one block for each sub-domain), should I? I tested the null space against the matrix using MatNullSpaceTest and the answer is true. Can I do something more to have a symmetric solution as outcome of the solver? Thank you in advance for any comments and hints. Best regards, Marco Cisternino

Hello Barry. I tried to force the solver to start from an initial guess which is not the solution of the problem. For sake of completeness, the solution has to be a constant field. With this initial condition, the solver iterates to a solution which is constant in the 2 sub-domains but * the constants have not the same value * they are not close to zero (minimal norm solution) * they are not opposite (zero-average solution over the whole domain, like 3 and -3) After 20 CFD iterations my pressure is 32 in one sub-domain and 2.2 in the other one. And their norm is increasing. How can I force the solver to give me minimal norm solution, or in other words the zero constant? I can do it by myself, anchoring domain-by-domain the solution removing its local average, but I was wondering if the solver can do this for me. In some way, giving a null space made of 2 vectors (1 on dofs living in the sub-domain and zero elsewhere), I would expect a solution with zero average in the 2 sub-domains, separately, but I’m wrong, probably. Finally, which is the closure of the problem defining the value of the constant? Zero-average condition, minimal norm condition, or none of them? Thanks! Bests, Marco Cisternino, PhD [email protected]<mailto:[email protected]> ______________________ Optimad Engineering Srl Via Bligny 5, Torino, Italia. +3901119719782 www.optimad.it<http://www.optimad.it/> From: Barry Smith <[email protected]> Sent: venerdì 1 ottobre 2021 16:56 To: Marco Cisternino <[email protected]> Cc: [email protected] Subject: Re: [petsc-users] Disconnected domains and Poisson equation On Oct 1, 2021, at 6:38 AM, Marco Cisternino <[email protected]<mailto:[email protected]>> wrote: Thank you Barry. I added a custom atoll = 1.0e-12 and this makes the CFD stable with all the linear solver types. CFD solution is good and pressure is a good “zero” field at every CFD iteration. I did the same test using ASM+ILU+FGMRES(BCGS and GMRES) and the behaviour is the same. During some CFD iteration the residual of linear system starts slightly higher than atol and the linear solver makes some iteration (2/3 iterations) before it stops because of atol. The pressure is still different in the 2 sub-domains (order 1.0e-14 because of those few linear solver iterations), therefore no symmetry of the solution In the 2 sub-domains. I think it is a matter of round-off, do you agree on this? Or do I need to take care of this difference as a symptom of something wrong? Yes, if the differences in the two solutions are order 1.e-14 that is very good, one cannot expect them to be identical. Thank you for your support. Marco Cisternino From: Barry Smith <[email protected]<mailto:[email protected]>> Sent: giovedì 30 settembre 2021 16:39 To: Marco Cisternino <[email protected]<mailto:[email protected]>> Cc: [email protected]<mailto:[email protected]> Subject: Re: [petsc-users] Disconnected domains and Poisson equation It looks like the initial solution (guess) is to round-off the solution to the linear system 9.010260489109e-14 0 KSP unpreconditioned resid norm 9.010260489109e-14 true resid norm 9.010260489109e-14 ||r(i)||/||b|| 2.021559024868e+00 0 KSP Residual norm 9.010260489109e-14 % max 1.000000000000e+00 min 1.000000000000e+00 max/min 1.000000000000e+00 1 KSP unpreconditioned resid norm 4.918108339808e-15 true resid norm 4.918171792537e-15 ||r(i)||/||b|| 1.103450292594e-01 1 KSP Residual norm 4.918108339808e-15 % max 9.566256813737e-01 min 9.566256813737e-01 max/min 1.000000000000e+00 2 KSP unpreconditioned resid norm 1.443599554690e-15 true resid norm 1.444867143493e-15 ||r(i)||/||b|| 3.241731154382e-02 2 KSP Residual norm 1.443599554690e-15 % max 9.614019380614e-01 min 7.360950481750e-01 max/min 1.306083963538e+00 Thus the Krylov solver will not be able to improve the solution, it then gets stuck trying to improve the solution but cannot because of round off. In other words the algorithm has converged (even at the initial solution (guess) and should stop immediately. You can use -ksp_atol 1.e-12 to get it to stop immediately without iterating if the initial residual is less than 1e-12. Barry On Sep 30, 2021, at 4:16 AM, Marco Cisternino <[email protected]<mailto:[email protected]>> wrote: Hello Barry. This is the output of ksp_view using fgmres and gamg. It has to be said that the solution of the linear system should be a zero values field. As you can see both unpreconditioned residual and r/b converge at this iteration of the CFD solver. During the time integration of the CFD, I can observe pressure linear solver residuals behaving in a different way: unpreconditioned residual stil converges but r/b stalls. After the output of ksp_view I add the output of ksp_monitor_true_residual for one of these iteration where r/b stalls. Thanks, KSP Object: 1 MPI processes type: fgmres restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement happy breakdown tolerance 1e-30 maximum iterations=100, nonzero initial guess tolerances: relative=1e-05, absolute=1e-50, divergence=10000. right preconditioning using UNPRECONDITIONED norm type for convergence test PC Object: 1 MPI processes type: gamg type is MULTIPLICATIVE, levels=4 cycles=v Cycles per PCApply=1 Using externally compute Galerkin coarse grid matrices GAMG specific options Threshold for dropping small values in graph on each level = 0.02 0.02 Threshold scaling factor for each level not specified = 1. AGG specific options Symmetric graph true Number of levels to square graph 1 Number smoothing steps 0 Coarse grid solver -- level ------------------------------- KSP Object: (mg_coarse_) 1 MPI processes type: preonly maximum iterations=10000, initial guess is zero tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using NONE norm type for convergence test PC Object: (mg_coarse_) 1 MPI processes type: bjacobi number of blocks = 1 Local solve is same for all blocks, in the following KSP and PC objects: KSP Object: (mg_coarse_sub_) 1 MPI processes type: preonly maximum iterations=1, initial guess is zero tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using DEFAULT norm type for convergence test PC Object: (mg_coarse_sub_) 1 MPI processes type: lu PC has not been set up so information may be incomplete out-of-place factorization tolerance for zero pivot 2.22045e-14 using diagonal shift on blocks to prevent zero pivot [INBLOCKS] matrix ordering: nd linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=18, cols=18 total: nonzeros=104, allocated nonzeros=104 total number of mallocs used during MatSetValues calls =0 not using I-node routines linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=18, cols=18 total: nonzeros=104, allocated nonzeros=104 total number of mallocs used during MatSetValues calls =0 not using I-node routines Down solver (pre-smoother) on level 1 ------------------------------- KSP Object: (mg_levels_1_) 1 MPI processes type: chebyshev eigenvalue estimates used: min = 0., max = 0. eigenvalues estimate via gmres min 0., max 0. eigenvalues estimated using gmres with translations [0. 0.1; 0. 1.1] KSP Object: (mg_levels_1_esteig_) 1 MPI processes type: gmres restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement happy breakdown tolerance 1e-30 maximum iterations=10, initial guess is zero tolerances: relative=1e-12, absolute=1e-50, divergence=10000. left preconditioning using DEFAULT norm type for convergence test estimating eigenvalues using noisy right hand side maximum iterations=2, nonzero initial guess tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using NONE norm type for convergence test PC Object: (mg_levels_1_) 1 MPI processes type: sor type = local_symmetric, iterations = 1, local iterations = 1, omega = 1. linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=67, cols=67 total: nonzeros=675, allocated nonzeros=675 total number of mallocs used during MatSetValues calls =0 not using I-node routines Up solver (post-smoother) same as down solver (pre-smoother) Down solver (pre-smoother) on level 2 ------------------------------- KSP Object: (mg_levels_2_) 1 MPI processes type: chebyshev eigenvalue estimates used: min = 0., max = 0. eigenvalues estimate via gmres min 0., max 0. eigenvalues estimated using gmres with translations [0. 0.1; 0. 1.1] KSP Object: (mg_levels_2_esteig_) 1 MPI processes type: gmres restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement happy breakdown tolerance 1e-30 maximum iterations=10, initial guess is zero tolerances: relative=1e-12, absolute=1e-50, divergence=10000. left preconditioning using DEFAULT norm type for convergence test estimating eigenvalues using noisy right hand side maximum iterations=2, nonzero initial guess tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using NONE norm type for convergence test PC Object: (mg_levels_2_) 1 MPI processes type: sor type = local_symmetric, iterations = 1, local iterations = 1, omega = 1. linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=348, cols=348 total: nonzeros=3928, allocated nonzeros=3928 total number of mallocs used during MatSetValues calls =0 not using I-node routines Up solver (post-smoother) same as down solver (pre-smoother) Down solver (pre-smoother) on level 3 ------------------------------- KSP Object: (mg_levels_3_) 1 MPI processes type: chebyshev eigenvalue estimates used: min = 0., max = 0. eigenvalues estimate via gmres min 0., max 0. eigenvalues estimated using gmres with translations [0. 0.1; 0. 1.1] KSP Object: (mg_levels_3_esteig_) 1 MPI processes type: gmres restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement happy breakdown tolerance 1e-30 maximum iterations=10, initial guess is zero tolerances: relative=1e-12, absolute=1e-50, divergence=10000. left preconditioning using DEFAULT norm type for convergence test estimating eigenvalues using noisy right hand side maximum iterations=2, nonzero initial guess tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using NONE norm type for convergence test PC Object: (mg_levels_3_) 1 MPI processes type: sor type = local_symmetric, iterations = 1, local iterations = 1, omega = 1. linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=3584, cols=3584 total: nonzeros=23616, allocated nonzeros=23616 total number of mallocs used during MatSetValues calls =0 has attached null space not using I-node routines Up solver (post-smoother) same as down solver (pre-smoother) linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=3584, cols=3584 total: nonzeros=23616, allocated nonzeros=23616 total number of mallocs used during MatSetValues calls =0 has attached null space not using I-node routines Pressure system has reached convergence in 0 iterations with reason 3. 0 KSP unpreconditioned resid norm 4.798763170703e-16 true resid norm 4.798763170703e-16 ||r(i)||/||b|| 1.000000000000e+00 0 KSP Residual norm 4.798763170703e-16 % max 1.000000000000e+00 min 1.000000000000e+00 max/min 1.000000000000e+00 1 KSP unpreconditioned resid norm 1.648749109132e-17 true resid norm 1.648749109132e-17 ||r(i)||/||b|| 3.435779284125e-02 1 KSP Residual norm 1.648749109132e-17 % max 9.561792537103e-01 min 9.561792537103e-01 max/min 1.000000000000e+00 2 KSP unpreconditioned resid norm 4.737880600040e-19 true resid norm 4.737880600040e-19 ||r(i)||/||b|| 9.873128619820e-04 2 KSP Residual norm 4.737880600040e-19 % max 9.828636644296e-01 min 9.293131521763e-01 max/min 1.057623753767e+00 3 KSP unpreconditioned resid norm 2.542212716830e-20 true resid norm 2.542212716830e-20 ||r(i)||/||b|| 5.297641551371e-05 3 KSP Residual norm 2.542212716830e-20 % max 9.933572357920e-01 min 9.158303248850e-01 max/min 1.084652046127e+00 4 KSP unpreconditioned resid norm 6.614510286263e-21 true resid norm 6.614510286269e-21 ||r(i)||/||b|| 1.378378146822e-05 4 KSP Residual norm 6.614510286263e-21 % max 9.950912550705e-01 min 6.296575800237e-01 max/min 1.580368896747e+00 5 KSP unpreconditioned resid norm 1.981505525281e-22 true resid norm 1.981505525272e-22 ||r(i)||/||b|| 4.129200493513e-07 5 KSP Residual norm 1.981505525281e-22 % max 9.984097962703e-01 min 5.316259535293e-01 max/min 1.878030577029e+00 Linear solve converged due to CONVERGED_RTOL iterations 5 Ksp_monitor_true_residual output for stalling r/b CFD iteration 0 KSP unpreconditioned resid norm 9.010260489109e-14 true resid norm 9.010260489109e-14 ||r(i)||/||b|| 2.021559024868e+00 0 KSP Residual norm 9.010260489109e-14 % max 1.000000000000e+00 min 1.000000000000e+00 max/min 1.000000000000e+00 1 KSP unpreconditioned resid norm 4.918108339808e-15 true resid norm 4.918171792537e-15 ||r(i)||/||b|| 1.103450292594e-01 1 KSP Residual norm 4.918108339808e-15 % max 9.566256813737e-01 min 9.566256813737e-01 max/min 1.000000000000e+00 2 KSP unpreconditioned resid norm 1.443599554690e-15 true resid norm 1.444867143493e-15 ||r(i)||/||b|| 3.241731154382e-02 2 KSP Residual norm 1.443599554690e-15 % max 9.614019380614e-01 min 7.360950481750e-01 max/min 1.306083963538e+00 3 KSP unpreconditioned resid norm 6.623206616803e-16 true resid norm 6.654132553541e-16 ||r(i)||/||b|| 1.492933720678e-02 3 KSP Residual norm 6.623206616803e-16 % max 9.764112945239e-01 min 4.911485418014e-01 max/min 1.988016274960e+00 4 KSP unpreconditioned resid norm 6.551896936698e-16 true resid norm 6.646157296305e-16 ||r(i)||/||b|| 1.491144376933e-02 4 KSP Residual norm 6.551896936698e-16 % max 9.883425885532e-01 min 1.461270778833e-01 max/min 6.763582786091e+00 5 KSP unpreconditioned resid norm 6.222297644887e-16 true resid norm 1.720560536914e-15 ||r(i)||/||b|| 3.860282047823e-02 5 KSP Residual norm 6.222297644887e-16 % max 1.000409371755e+00 min 4.989767363560e-03 max/min 2.004921870829e+02 6 KSP unpreconditioned resid norm 6.496945794974e-17 true resid norm 2.031914800253e-14 ||r(i)||/||b|| 4.558842341106e-01 6 KSP Residual norm 6.496945794974e-17 % max 1.004914985753e+00 min 1.459258738706e-03 max/min 6.886475709192e+02 7 KSP unpreconditioned resid norm 1.965237342540e-17 true resid norm 1.684522207337e-14 ||r(i)||/||b|| 3.779425772373e-01 7 KSP Residual norm 1.965237342540e-17 % max 1.005737762541e+00 min 1.452603803766e-03 max/min 6.923689446035e+02 8 KSP unpreconditioned resid norm 1.627718951285e-17 true resid norm 1.958642967520e-14 ||r(i)||/||b|| 4.394448276241e-01 8 KSP Residual norm 1.627718951285e-17 % max 1.006364278765e+00 min 1.452081813014e-03 max/min 6.930492963590e+02 9 KSP unpreconditioned resid norm 1.616577677764e-17 true resid norm 2.019110946644e-14 ||r(i)||/||b|| 4.530115373837e-01 9 KSP Residual norm 1.616577677764e-17 % max 1.006648747131e+00 min 1.452031376577e-03 max/min 6.932692801059e+02 10 KSP unpreconditioned resid norm 1.285788988203e-17 true resid norm 2.065082694477e-14 ||r(i)||/||b|| 4.633258453698e-01 10 KSP Residual norm 1.285788988203e-17 % max 1.007469033514e+00 min 1.433291867068e-03 max/min 7.029057072477e+02 11 KSP unpreconditioned resid norm 5.490854431580e-19 true resid norm 1.798071628891e-14 ||r(i)||/||b|| 4.034187394623e-01 11 KSP Residual norm 5.490854431580e-19 % max 1.008058905554e+00 min 1.369401685301e-03 max/min 7.361309076612e+02 12 KSP unpreconditioned resid norm 1.371754802104e-20 true resid norm 1.965688920064e-14 ||r(i)||/||b|| 4.410256708163e-01 12 KSP Residual norm 1.371754802104e-20 % max 1.008409402214e+00 min 1.369243011779e-03 max/min 7.364721919624e+02 Linear solve converged due to CONVERGED_RTOL iterations 12 Marco Cisternino From: Barry Smith <[email protected]<mailto:[email protected]>> Sent: mercoledì 29 settembre 2021 18:34 To: Marco Cisternino <[email protected]<mailto:[email protected]>> Cc: [email protected]<mailto:[email protected]> Subject: Re: [petsc-users] Disconnected domains and Poisson equation On Sep 29, 2021, at 11:59 AM, Marco Cisternino <[email protected]<mailto:[email protected]>> wrote: For sake of completeness, explicitly building the null space using a vector per sub-domain make s the CFD runs using BCGS and GMRES more stable, but still slower than FGMRES. Something is strange. Please run with -ksp_view and send the output on the solver details. I had divergence using BCGS and GMRES setting the null space with only one constant. Thanks Marco Cisternino From: Marco Cisternino Sent: mercoledì 29 settembre 2021 17:54 To: Barry Smith <[email protected]<mailto:[email protected]>> Cc: [email protected]<mailto:[email protected]> Subject: RE: [petsc-users] Disconnected domains and Poisson equation Thank you Barry for the quick reply. About the null space: I already tried what you suggest, building 2 Vec (constants) with 0 and 1 chosen by sub-domain, normalizing them and setting the null space like this MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,nconstants,constants,&nullspace); The solution is slightly different in values but it is still different in the two sub-domains. About the solver: I tried BCGS, GMRES and FGMRES. The linear system is a pressure system in a navier-stokes solver and only solving with FGMRES makes the CFD stable, with BCGS and GMRES the CFD solution diverges. Moreover, in the same case but with a single domain, CFD solution is stable using all the solvers, but FGMRES converges in much less iterations than the others. Marco Cisternino From: Barry Smith <[email protected]<mailto:[email protected]>> Sent: mercoledì 29 settembre 2021 15:59 To: Marco Cisternino <[email protected]<mailto:[email protected]>> Cc: [email protected]<mailto:[email protected]> Subject: Re: [petsc-users] Disconnected domains and Poisson equation The problem actually has a two dimensional null space; constant on each domain but possibly different constants. I think you need to build the MatNullSpace by explicitly constructing two vectors, one with 0 on one domain and constant value on the other and one with 0 on the other domain and constant on the first. Separate note: why use FGMRES instead of just GMRES? If the problem is linear and the preconditioner is linear (no GMRES inside the smoother) then you can just use GMRES and it will save a little space/work and be conceptually clearer. Barry On Sep 29, 2021, at 8:46 AM, Marco Cisternino <[email protected]<mailto:[email protected]>> wrote: Good morning, I want to solve the Poisson equation on a 3D domain with 2 non-connected sub-domains. I am using FGMRES+GAMG and I have no problem if the two sub-domains see a Dirichlet boundary condition each. On the same domain I would like to solve the Poisson equation imposing periodic boundary condition in one direction and homogenous Neumann boundary conditions in the other two directions. The two sub-domains are symmetric with respect to the separation between them and the operator discretization and the right hand side are symmetric as well. It would be nice to have the same solution in both the sub-domains. Setting the null space to the constant, the solver converges to a solution having the same gradients in both sub-domains but different values. Am I doing some wrong with the null space? I’m not setting a block matrix (one block for each sub-domain), should I? I tested the null space against the matrix using MatNullSpaceTest and the answer is true. Can I do something more to have a symmetric solution as outcome of the solver? Thank you in advance for any comments and hints. Best regards, Marco Cisternino

On Oct 6, 2021, at 1:20 PM, Marco Cisternino <[email protected]> wrote:

Hello Barry. I tried to force the solver to start from an initial guess which is not the solution of the problem. For sake of completeness, the solution has to be a constant field. With this initial condition, the solver iterates to a solution which is constant in the 2 sub-domains but the constants have not the same value they are not close to zero (minimal norm solution) they are not opposite (zero-average solution over the whole domain, like 3 and -3) After 20 CFD iterations my pressure is 32 in one sub-domain and 2.2 in the other one. And their norm is increasing. How can I force the solver to give me minimal norm solution, or in other words the zero constant?

Providing the appropriate null space should result in GMRES giving you the minimal norm solution which corresponds to the average of the solution on each domain being zero. For a general right hand side the residual will not decrease to zero because the right hand side is inconsistent. You can still run GMRES and the solution will converge but you will have a nonzero residual at the end (this makes the stopping criteria harder so it is useful to remove the inconsistent part of the right hand side out from the right hand side before calling GMRES.) Barry

I can do it by myself, anchoring domain-by-domain the solution removing its local average, but I was wondering if the solver can do this for me. In some way, giving a null space made of 2 vectors (1 on dofs living in the sub-domain and zero elsewhere), I would expect a solution with zero average in the 2 sub-domains, separately, but I’m wrong, probably. Finally, which is the closure of the problem defining the value of the constant? Zero-average condition, minimal norm condition, or none of them?

Thanks!

Bests,

Marco Cisternino, PhD [email protected] <mailto:[email protected]> ______________________ Optimad Engineering Srl Via Bligny 5, Torino, Italia. +3901119719782 www.optimad.it <http://www.optimad.it/>

From: Barry Smith <[email protected] <mailto:[email protected]>> Sent: venerdì 1 ottobre 2021 16:56 To: Marco Cisternino <[email protected] <mailto:[email protected]>> Cc: [email protected] <mailto:[email protected]> Subject: Re: [petsc-users] Disconnected domains and Poisson equation

On Oct 1, 2021, at 6:38 AM, Marco Cisternino <[email protected] <mailto:[email protected]>> wrote:

Thank you Barry. I added a custom atoll = 1.0e-12 and this makes the CFD stable with all the linear solver types. CFD solution is good and pressure is a good “zero” field at every CFD iteration. I did the same test using ASM+ILU+FGMRES(BCGS and GMRES) and the behaviour is the same. During some CFD iteration the residual of linear system starts slightly higher than atol and the linear solver makes some iteration (2/3 iterations) before it stops because of atol. The pressure is still different in the 2 sub-domains (order 1.0e-14 because of those few linear solver iterations), therefore no symmetry of the solution In the 2 sub-domains. I think it is a matter of round-off, do you agree on this? Or do I need to take care of this difference as a symptom of something wrong?

Yes, if the differences in the two solutions are order 1.e-14 that is very good, one cannot expect them to be identical.

Thank you for your support.

Marco Cisternino

From: Barry Smith <[email protected] <mailto:[email protected]>> Sent: giovedì 30 settembre 2021 16:39 To: Marco Cisternino <[email protected] <mailto:[email protected]>> Cc: [email protected] <mailto:[email protected]> Subject: Re: [petsc-users] Disconnected domains and Poisson equation

It looks like the initial solution (guess) is to round-off the solution to the linear system 9.010260489109e-14

0 KSP unpreconditioned resid norm 9.010260489109e-14 true resid norm 9.010260489109e-14 ||r(i)||/||b|| 2.021559024868e+00 0 KSP Residual norm 9.010260489109e-14 % max 1.000000000000e+00 min 1.000000000000e+00 max/min 1.000000000000e+00 1 KSP unpreconditioned resid norm 4.918108339808e-15 true resid norm 4.918171792537e-15 ||r(i)||/||b|| 1.103450292594e-01 1 KSP Residual norm 4.918108339808e-15 % max 9.566256813737e-01 min 9.566256813737e-01 max/min 1.000000000000e+00 2 KSP unpreconditioned resid norm 1.443599554690e-15 true resid norm 1.444867143493e-15 ||r(i)||/||b|| 3.241731154382e-02 2 KSP Residual norm 1.443599554690e-15 % max 9.614019380614e-01 min 7.360950481750e-01 max/min 1.306083963538e+00

Thus the Krylov solver will not be able to improve the solution, it then gets stuck trying to improve the solution but cannot because of round off.

In other words the algorithm has converged (even at the initial solution (guess) and should stop immediately.

You can use -ksp_atol 1.e-12 to get it to stop immediately without iterating if the initial residual is less than 1e-12.

Barry

On Sep 30, 2021, at 4:16 AM, Marco Cisternino <[email protected] <mailto:[email protected]>> wrote:

Hello Barry. This is the output of ksp_view using fgmres and gamg. It has to be said that the solution of the linear system should be a zero values field. As you can see both unpreconditioned residual and r/b converge at this iteration of the CFD solver. During the time integration of the CFD, I can observe pressure linear solver residuals behaving in a different way: unpreconditioned residual stil converges but r/b stalls. After the output of ksp_view I add the output of ksp_monitor_true_residual for one of these iteration where r/b stalls. Thanks,

KSP Object: 1 MPI processes type: fgmres restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement happy breakdown tolerance 1e-30 maximum iterations=100, nonzero initial guess tolerances: relative=1e-05, absolute=1e-50, divergence=10000. right preconditioning using UNPRECONDITIONED norm type for convergence test PC Object: 1 MPI processes type: gamg type is MULTIPLICATIVE, levels=4 cycles=v Cycles per PCApply=1 Using externally compute Galerkin coarse grid matrices GAMG specific options Threshold for dropping small values in graph on each level = 0.02 0.02 Threshold scaling factor for each level not specified = 1. AGG specific options Symmetric graph true Number of levels to square graph 1 Number smoothing steps 0 Coarse grid solver -- level ------------------------------- KSP Object: (mg_coarse_) 1 MPI processes type: preonly maximum iterations=10000, initial guess is zero tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using NONE norm type for convergence test PC Object: (mg_coarse_) 1 MPI processes type: bjacobi number of blocks = 1 Local solve is same for all blocks, in the following KSP and PC objects: KSP Object: (mg_coarse_sub_) 1 MPI processes type: preonly maximum iterations=1, initial guess is zero tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using DEFAULT norm type for convergence test PC Object: (mg_coarse_sub_) 1 MPI processes type: lu PC has not been set up so information may be incomplete out-of-place factorization tolerance for zero pivot 2.22045e-14 using diagonal shift on blocks to prevent zero pivot [INBLOCKS] matrix ordering: nd linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=18, cols=18 total: nonzeros=104, allocated nonzeros=104 total number of mallocs used during MatSetValues calls =0 not using I-node routines linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=18, cols=18 total: nonzeros=104, allocated nonzeros=104 total number of mallocs used during MatSetValues calls =0 not using I-node routines Down solver (pre-smoother) on level 1 ------------------------------- KSP Object: (mg_levels_1_) 1 MPI processes type: chebyshev eigenvalue estimates used: min = 0., max = 0. eigenvalues estimate via gmres min 0., max 0. eigenvalues estimated using gmres with translations [0. 0.1; 0. 1.1] KSP Object: (mg_levels_1_esteig_) 1 MPI processes type: gmres restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement happy breakdown tolerance 1e-30 maximum iterations=10, initial guess is zero tolerances: relative=1e-12, absolute=1e-50, divergence=10000. left preconditioning using DEFAULT norm type for convergence test estimating eigenvalues using noisy right hand side maximum iterations=2, nonzero initial guess tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using NONE norm type for convergence test PC Object: (mg_levels_1_) 1 MPI processes type: sor type = local_symmetric, iterations = 1, local iterations = 1, omega = 1. linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=67, cols=67 total: nonzeros=675, allocated nonzeros=675 total number of mallocs used during MatSetValues calls =0 not using I-node routines Up solver (post-smoother) same as down solver (pre-smoother) Down solver (pre-smoother) on level 2 ------------------------------- KSP Object: (mg_levels_2_) 1 MPI processes type: chebyshev eigenvalue estimates used: min = 0., max = 0. eigenvalues estimate via gmres min 0., max 0. eigenvalues estimated using gmres with translations [0. 0.1; 0. 1.1] KSP Object: (mg_levels_2_esteig_) 1 MPI processes type: gmres restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement happy breakdown tolerance 1e-30 maximum iterations=10, initial guess is zero tolerances: relative=1e-12, absolute=1e-50, divergence=10000. left preconditioning using DEFAULT norm type for convergence test estimating eigenvalues using noisy right hand side maximum iterations=2, nonzero initial guess tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using NONE norm type for convergence test PC Object: (mg_levels_2_) 1 MPI processes type: sor type = local_symmetric, iterations = 1, local iterations = 1, omega = 1. linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=348, cols=348 total: nonzeros=3928, allocated nonzeros=3928 total number of mallocs used during MatSetValues calls =0 not using I-node routines Up solver (post-smoother) same as down solver (pre-smoother) Down solver (pre-smoother) on level 3 ------------------------------- KSP Object: (mg_levels_3_) 1 MPI processes type: chebyshev eigenvalue estimates used: min = 0., max = 0. eigenvalues estimate via gmres min 0., max 0. eigenvalues estimated using gmres with translations [0. 0.1; 0. 1.1] KSP Object: (mg_levels_3_esteig_) 1 MPI processes type: gmres restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement happy breakdown tolerance 1e-30 maximum iterations=10, initial guess is zero tolerances: relative=1e-12, absolute=1e-50, divergence=10000. left preconditioning using DEFAULT norm type for convergence test estimating eigenvalues using noisy right hand side maximum iterations=2, nonzero initial guess tolerances: relative=1e-05, absolute=1e-50, divergence=10000. left preconditioning using NONE norm type for convergence test PC Object: (mg_levels_3_) 1 MPI processes type: sor type = local_symmetric, iterations = 1, local iterations = 1, omega = 1. linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=3584, cols=3584 total: nonzeros=23616, allocated nonzeros=23616 total number of mallocs used during MatSetValues calls =0 has attached null space not using I-node routines Up solver (post-smoother) same as down solver (pre-smoother) linear system matrix = precond matrix: Mat Object: 1 MPI processes type: seqaij rows=3584, cols=3584 total: nonzeros=23616, allocated nonzeros=23616 total number of mallocs used during MatSetValues calls =0 has attached null space not using I-node routines Pressure system has reached convergence in 0 iterations with reason 3. 0 KSP unpreconditioned resid norm 4.798763170703e-16 true resid norm 4.798763170703e-16 ||r(i)||/||b|| 1.000000000000e+00 0 KSP Residual norm 4.798763170703e-16 % max 1.000000000000e+00 min 1.000000000000e+00 max/min 1.000000000000e+00 1 KSP unpreconditioned resid norm 1.648749109132e-17 true resid norm 1.648749109132e-17 ||r(i)||/||b|| 3.435779284125e-02 1 KSP Residual norm 1.648749109132e-17 % max 9.561792537103e-01 min 9.561792537103e-01 max/min 1.000000000000e+00 2 KSP unpreconditioned resid norm 4.737880600040e-19 true resid norm 4.737880600040e-19 ||r(i)||/||b|| 9.873128619820e-04 2 KSP Residual norm 4.737880600040e-19 % max 9.828636644296e-01 min 9.293131521763e-01 max/min 1.057623753767e+00 3 KSP unpreconditioned resid norm 2.542212716830e-20 true resid norm 2.542212716830e-20 ||r(i)||/||b|| 5.297641551371e-05 3 KSP Residual norm 2.542212716830e-20 % max 9.933572357920e-01 min 9.158303248850e-01 max/min 1.084652046127e+00 4 KSP unpreconditioned resid norm 6.614510286263e-21 true resid norm 6.614510286269e-21 ||r(i)||/||b|| 1.378378146822e-05 4 KSP Residual norm 6.614510286263e-21 % max 9.950912550705e-01 min 6.296575800237e-01 max/min 1.580368896747e+00 5 KSP unpreconditioned resid norm 1.981505525281e-22 true resid norm 1.981505525272e-22 ||r(i)||/||b|| 4.129200493513e-07 5 KSP Residual norm 1.981505525281e-22 % max 9.984097962703e-01 min 5.316259535293e-01 max/min 1.878030577029e+00 Linear solve converged due to CONVERGED_RTOL iterations 5

Ksp_monitor_true_residual output for stalling r/b CFD iteration 0 KSP unpreconditioned resid norm 9.010260489109e-14 true resid norm 9.010260489109e-14 ||r(i)||/||b|| 2.021559024868e+00 0 KSP Residual norm 9.010260489109e-14 % max 1.000000000000e+00 min 1.000000000000e+00 max/min 1.000000000000e+00 1 KSP unpreconditioned resid norm 4.918108339808e-15 true resid norm 4.918171792537e-15 ||r(i)||/||b|| 1.103450292594e-01 1 KSP Residual norm 4.918108339808e-15 % max 9.566256813737e-01 min 9.566256813737e-01 max/min 1.000000000000e+00 2 KSP unpreconditioned resid norm 1.443599554690e-15 true resid norm 1.444867143493e-15 ||r(i)||/||b|| 3.241731154382e-02 2 KSP Residual norm 1.443599554690e-15 % max 9.614019380614e-01 min 7.360950481750e-01 max/min 1.306083963538e+00 3 KSP unpreconditioned resid norm 6.623206616803e-16 true resid norm 6.654132553541e-16 ||r(i)||/||b|| 1.492933720678e-02 3 KSP Residual norm 6.623206616803e-16 % max 9.764112945239e-01 min 4.911485418014e-01 max/min 1.988016274960e+00 4 KSP unpreconditioned resid norm 6.551896936698e-16 true resid norm 6.646157296305e-16 ||r(i)||/||b|| 1.491144376933e-02 4 KSP Residual norm 6.551896936698e-16 % max 9.883425885532e-01 min 1.461270778833e-01 max/min 6.763582786091e+00 5 KSP unpreconditioned resid norm 6.222297644887e-16 true resid norm 1.720560536914e-15 ||r(i)||/||b|| 3.860282047823e-02 5 KSP Residual norm 6.222297644887e-16 % max 1.000409371755e+00 min 4.989767363560e-03 max/min 2.004921870829e+02 6 KSP unpreconditioned resid norm 6.496945794974e-17 true resid norm 2.031914800253e-14 ||r(i)||/||b|| 4.558842341106e-01 6 KSP Residual norm 6.496945794974e-17 % max 1.004914985753e+00 min 1.459258738706e-03 max/min 6.886475709192e+02 7 KSP unpreconditioned resid norm 1.965237342540e-17 true resid norm 1.684522207337e-14 ||r(i)||/||b|| 3.779425772373e-01 7 KSP Residual norm 1.965237342540e-17 % max 1.005737762541e+00 min 1.452603803766e-03 max/min 6.923689446035e+02 8 KSP unpreconditioned resid norm 1.627718951285e-17 true resid norm 1.958642967520e-14 ||r(i)||/||b|| 4.394448276241e-01 8 KSP Residual norm 1.627718951285e-17 % max 1.006364278765e+00 min 1.452081813014e-03 max/min 6.930492963590e+02 9 KSP unpreconditioned resid norm 1.616577677764e-17 true resid norm 2.019110946644e-14 ||r(i)||/||b|| 4.530115373837e-01 9 KSP Residual norm 1.616577677764e-17 % max 1.006648747131e+00 min 1.452031376577e-03 max/min 6.932692801059e+02 10 KSP unpreconditioned resid norm 1.285788988203e-17 true resid norm 2.065082694477e-14 ||r(i)||/||b|| 4.633258453698e-01 10 KSP Residual norm 1.285788988203e-17 % max 1.007469033514e+00 min 1.433291867068e-03 max/min 7.029057072477e+02 11 KSP unpreconditioned resid norm 5.490854431580e-19 true resid norm 1.798071628891e-14 ||r(i)||/||b|| 4.034187394623e-01 11 KSP Residual norm 5.490854431580e-19 % max 1.008058905554e+00 min 1.369401685301e-03 max/min 7.361309076612e+02 12 KSP unpreconditioned resid norm 1.371754802104e-20 true resid norm 1.965688920064e-14 ||r(i)||/||b|| 4.410256708163e-01 12 KSP Residual norm 1.371754802104e-20 % max 1.008409402214e+00 min 1.369243011779e-03 max/min 7.364721919624e+02 Linear solve converged due to CONVERGED_RTOL iterations 12

Marco Cisternino

From: Barry Smith <[email protected] <mailto:[email protected]>> Sent: mercoledì 29 settembre 2021 18:34 To: Marco Cisternino <[email protected] <mailto:[email protected]>> Cc: [email protected] <mailto:[email protected]> Subject: Re: [petsc-users] Disconnected domains and Poisson equation

On Sep 29, 2021, at 11:59 AM, Marco Cisternino <[email protected] <mailto:[email protected]>> wrote:

For sake of completeness, explicitly building the null space using a vector per sub-domain make s the CFD runs using BCGS and GMRES more stable, but still slower than FGMRES.

Something is strange. Please run with -ksp_view and send the output on the solver details.

I had divergence using BCGS and GMRES setting the null space with only one constant. Thanks

Marco Cisternino

From: Marco Cisternino Sent: mercoledì 29 settembre 2021 17:54 To: Barry Smith <[email protected] <mailto:[email protected]>> Cc: [email protected] <mailto:[email protected]> Subject: RE: [petsc-users] Disconnected domains and Poisson equation

Thank you Barry for the quick reply. About the null space: I already tried what you suggest, building 2 Vec (constants) with 0 and 1 chosen by sub-domain, normalizing them and setting the null space like this MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,nconstants,constants,&nullspace); The solution is slightly different in values but it is still different in the two sub-domains. About the solver: I tried BCGS, GMRES and FGMRES. The linear system is a pressure system in a navier-stokes solver and only solving with FGMRES makes the CFD stable, with BCGS and GMRES the CFD solution diverges. Moreover, in the same case but with a single domain, CFD solution is stable using all the solvers, but FGMRES converges in much less iterations than the others.

Marco Cisternino

From: Barry Smith <[email protected] <mailto:[email protected]>> Sent: mercoledì 29 settembre 2021 15:59 To: Marco Cisternino <[email protected] <mailto:[email protected]>> Cc: [email protected] <mailto:[email protected]> Subject: Re: [petsc-users] Disconnected domains and Poisson equation

The problem actually has a two dimensional null space; constant on each domain but possibly different constants. I think you need to build the MatNullSpace by explicitly constructing two vectors, one with 0 on one domain and constant value on the other and one with 0 on the other domain and constant on the first.

Separate note: why use FGMRES instead of just GMRES? If the problem is linear and the preconditioner is linear (no GMRES inside the smoother) then you can just use GMRES and it will save a little space/work and be conceptually clearer.

Barry

On Sep 29, 2021, at 8:46 AM, Marco Cisternino <[email protected] <mailto:[email protected]>> wrote:

Good morning, I want to solve the Poisson equation on a 3D domain with 2 non-connected sub-domains. I am using FGMRES+GAMG and I have no problem if the two sub-domains see a Dirichlet boundary condition each. On the same domain I would like to solve the Poisson equation imposing periodic boundary condition in one direction and homogenous Neumann boundary conditions in the other two directions. The two sub-domains are symmetric with respect to the separation between them and the operator discretization and the right hand side are symmetric as well. It would be nice to have the same solution in both the sub-domains. Setting the null space to the constant, the solver converges to a solution having the same gradients in both sub-domains but different values. Am I doing some wrong with the null space? I’m not setting a block matrix (one block for each sub-domain), should I? I tested the null space against the matrix using MatNullSpaceTest and the answer is true. Can I do something more to have a symmetric solution as outcome of the solver? Thank you in advance for any comments and hints.

Best regards,

Marco Cisternino