Re: [petsc-users] finite difference jacobian errors when given non-constant initial condition

[Zou, Ling via petsc-users]

The other symptom is the same:

  *   Using coloring, finite differencing respects the specified non-zero 
pattern, but gives wrong (very large) Jacobian entries (J_ij)
  *   Using dense matrix assumption, finite differencing does not respect the 
non-zero pattern determined by your numeric, which is a clear sign of residual 
function code bug (your residual function does not respect your numeric).
-Ling

From: petsc-users  on behalf of Zou, Ling via 
petsc-users 
Date: Sunday, April 21, 2024 at 2:28 PM
To: Mark Lohry 
Cc: PETSc 
Subject: Re: [petsc-users] finite difference jacobian errors when given 
non-constant initial condition
Very interesting.
I happened to encounter something very similar a couple of days ago, which, of 
course, was due to a code bug I introduced. The code bug was in the residual 
function. I used a local vector to track ‘heat flux’, which should be zero-ed 
out at the beginning of each residual function evaluation. I did not zero it, 
and I observed very similar results, the Jacobian is completely wrong, with 
large values (J_ij keeps increasing after each iteration), and non-zero values 
are observed in locations which should be perfect zero. The symptom is very 
much like what you are seeing here. I suspect a similar bug. (Maybe you forgot 
zero the coefficients of P1 re-construction? Using constant value 1, 
reconstructed dphi/dx = 0, so however many iterations, still zero).

-Ling

From: Mark Lohry 
Date: Sunday, April 21, 2024 at 12:35 PM
To: Zou, Ling 
Cc: PETSc 
Subject: Re: [petsc-users] finite difference jacobian errors when given 
non-constant initial condition
The coloring I'm fairly confident is correct -- I use the same process for 3D 
unstructured grids and everything seems to work. The residual function is also 
validated. As a test I did as you suggested -- assume the matrix is dense -- and
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd
The coloring I'm fairly confident is correct -- I use the same process for 3D 
unstructured grids and everything seems to work. The residual function is also 
validated.

As a test I did as you suggested -- assume the matrix is dense -- and I get the 
same bad results, just now the zero blocks are filled.

Assuming dense, giving it a constant vector, all is good:

 4.23516e-16 -1.10266  0.31831   -0.085290900   
  -0.31831  1.18795
 1.10266 -4.23516e-16 -1.18795  0.3183100   
 0.0852909 -0.31831
-0.31831  1.18795  2.11758e-16 -1.10266  0.31831   -0.0852909   
 00
   0.0852909 -0.31831  1.10266 -4.23516e-16 -1.18795  0.31831   
 00
   00 -0.31831  1.18795  2.11758e-16 -1.10266   
   0.31831   -0.0852909
   000.0852909 -0.31831  1.10266 -4.23516e-16   
  -1.18795  0.31831
 0.31831   -0.085290900 -0.31831  1.18795  
4.23516e-16 -1.10266


Assuming dense, giving it sin(x), all is bad:

-1.76177e+08 -6.07287e+07 -6.07287e+07 -1.76177e+08  1.76177e+08  6.07287e+07  
6.07287e+07  1.76177e+08
-1.31161e+08 -4.52116e+07 -4.52116e+07 -1.31161e+08  1.31161e+08  4.52116e+07  
4.52116e+07  1.31161e+08
 1.31161e+08  4.52116e+07  4.52116e+07  1.31161e+08 -1.31161e+08 -4.52116e+07 
-4.52116e+07 -1.31161e+08
 1.76177e+08  6.07287e+07  6.07287e+07  1.76177e+08 -1.76177e+08 -6.07287e+07 
-6.07287e+07 -1.76177e+08
 1.76177e+08  6.07287e+07  6.07287e+07  1.76177e+08 -1.76177e+08 -6.07287e+07 
-6.07287e+07 -1.76177e+08
 1.31161e+08  4.52116e+07  4.52116e+07  1.31161e+08 -1.31161e+08 -4.52116e+07 
-4.52116e+07 -1.31161e+08
-1.31161e+08 -4.52116e+07 -4.52116e+07 -1.31161e+08  1.31161e+08  4.52116e+07  
4.52116e+07  1.31161e+08
-1.76177e+08 -6.07287e+07 -6.07287e+07 -1.76177e+08  1.76177e+08  6.07287e+07  
6.07287e+07  1.76177e+08

Scratching my head over here... I've been using these routines successfully for 
years in much more complex code.

On Sun, Apr 21, 2024 at 12:36 PM Zou, Ling mailto:l...@anl.gov>> 
wrote:
Edit:

  *   how do you do the coloring when using PETSc finite differencing? An 
incorrect coloring may give you wrong Jacobian. For debugging purpose, the 
simplest way to avoid an incorrect coloring is to assume the matrix is dense 
(slow but error proofing). If the numeric converges as expected, then fine tune 
your coloring to make it right and fast.


From: petsc-users 
mailto:petsc-users-boun...@mcs.anl.gov>> on 
behalf of Zou, Ling via petsc-users 
mailto:petsc-users@mcs.anl.gov>>
Date: Sunday, April 21, 2024 at 11:29 AM
To: Mark Lohry mailto:mlo...@gmail.com>>, PETSc 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] finite difference jacobian errors when given 
non-constant initial condition
Hi Mark, I am working on a project having similar numeric you

Re: [petsc-users] finite difference jacobian errors when given non-constant initial condition

[Zou, Ling via petsc-users]

Very interesting.
I happened to encounter something very similar a couple of days ago, which, of 
course, was due to a code bug I introduced. The code bug was in the residual 
function. I used a local vector to track ‘heat flux’, which should be zero-ed 
out at the beginning of each residual function evaluation. I did not zero it, 
and I observed very similar results, the Jacobian is completely wrong, with 
large values (J_ij keeps increasing after each iteration), and non-zero values 
are observed in locations which should be perfect zero. The symptom is very 
much like what you are seeing here. I suspect a similar bug. (Maybe you forgot 
zero the coefficients of P1 re-construction? Using constant value 1, 
reconstructed dphi/dx = 0, so however many iterations, still zero).

-Ling

From: Mark Lohry 
Date: Sunday, April 21, 2024 at 12:35 PM
To: Zou, Ling 
Cc: PETSc 
Subject: Re: [petsc-users] finite difference jacobian errors when given 
non-constant initial condition
The coloring I'm fairly confident is correct -- I use the same process for 3D 
unstructured grids and everything seems to work. The residual function is also 
validated. As a test I did as you suggested -- assume the matrix is dense -- and
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd
The coloring I'm fairly confident is correct -- I use the same process for 3D 
unstructured grids and everything seems to work. The residual function is also 
validated.

As a test I did as you suggested -- assume the matrix is dense -- and I get the 
same bad results, just now the zero blocks are filled.

Assuming dense, giving it a constant vector, all is good:

 4.23516e-16 -1.10266  0.31831   -0.085290900   
  -0.31831  1.18795
 1.10266 -4.23516e-16 -1.18795  0.3183100   
 0.0852909 -0.31831
-0.31831  1.18795  2.11758e-16 -1.10266  0.31831   -0.0852909   
 00
   0.0852909 -0.31831  1.10266 -4.23516e-16 -1.18795  0.31831   
 00
   00 -0.31831  1.18795  2.11758e-16 -1.10266   
   0.31831   -0.0852909
   000.0852909 -0.31831  1.10266 -4.23516e-16   
  -1.18795  0.31831
 0.31831   -0.085290900 -0.31831  1.18795  
4.23516e-16 -1.10266


Assuming dense, giving it sin(x), all is bad:

-1.76177e+08 -6.07287e+07 -6.07287e+07 -1.76177e+08  1.76177e+08  6.07287e+07  
6.07287e+07  1.76177e+08
-1.31161e+08 -4.52116e+07 -4.52116e+07 -1.31161e+08  1.31161e+08  4.52116e+07  
4.52116e+07  1.31161e+08
 1.31161e+08  4.52116e+07  4.52116e+07  1.31161e+08 -1.31161e+08 -4.52116e+07 
-4.52116e+07 -1.31161e+08
 1.76177e+08  6.07287e+07  6.07287e+07  1.76177e+08 -1.76177e+08 -6.07287e+07 
-6.07287e+07 -1.76177e+08
 1.76177e+08  6.07287e+07  6.07287e+07  1.76177e+08 -1.76177e+08 -6.07287e+07 
-6.07287e+07 -1.76177e+08
 1.31161e+08  4.52116e+07  4.52116e+07  1.31161e+08 -1.31161e+08 -4.52116e+07 
-4.52116e+07 -1.31161e+08
-1.31161e+08 -4.52116e+07 -4.52116e+07 -1.31161e+08  1.31161e+08  4.52116e+07  
4.52116e+07  1.31161e+08
-1.76177e+08 -6.07287e+07 -6.07287e+07 -1.76177e+08  1.76177e+08  6.07287e+07  
6.07287e+07  1.76177e+08

Scratching my head over here... I've been using these routines successfully for 
years in much more complex code.

On Sun, Apr 21, 2024 at 12:36 PM Zou, Ling mailto:l...@anl.gov>> 
wrote:
Edit:

  *   how do you do the coloring when using PETSc finite differencing? An 
incorrect coloring may give you wrong Jacobian. For debugging purpose, the 
simplest way to avoid an incorrect coloring is to assume the matrix is dense 
(slow but error proofing). If the numeric converges as expected, then fine tune 
your coloring to make it right and fast.


From: petsc-users 
mailto:petsc-users-boun...@mcs.anl.gov>> on 
behalf of Zou, Ling via petsc-users 
mailto:petsc-users@mcs.anl.gov>>
Date: Sunday, April 21, 2024 at 11:29 AM
To: Mark Lohry mailto:mlo...@gmail.com>>, PETSc 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] finite difference jacobian errors when given 
non-constant initial condition
Hi Mark, I am working on a project having similar numeric you have, 
one-dimensional finite volume method with second-order slope limiter TVD, and 
PETSc finite differencing gives perfect Jacobian even for complex problems.
So, I tend to believe that your implementation may have some problem. Some 
lessons I learned during my code development:


  *   how do you do the coloring when using PETSc finite differencing? An 
incorrect coloring may give you wrong Jacobian. The simplest way to avoid an 
incorrect coloring is to assume the matrix is dense (slow but error proofing).
  *   Residual function evaluation not correctly implemented can also lead to 
incorrect Jacobian. In your case, you m

Re: [petsc-users] finite difference jacobian errors when given non-constant initial condition

[Zou, Ling via petsc-users]

Edit:

  *   how do you do the coloring when using PETSc finite differencing? An 
incorrect coloring may give you wrong Jacobian. For debugging purpose, the 
simplest way to avoid an incorrect coloring is to assume the matrix is dense 
(slow but error proofing). If the numeric converges as expected, then fine tune 
your coloring to make it right and fast.


From: petsc-users  on behalf of Zou, Ling via 
petsc-users 
Date: Sunday, April 21, 2024 at 11:29 AM
To: Mark Lohry , PETSc 
Subject: Re: [petsc-users] finite difference jacobian errors when given 
non-constant initial condition
Hi Mark, I am working on a project having similar numeric you have, 
one-dimensional finite volume method with second-order slope limiter TVD, and 
PETSc finite differencing gives perfect Jacobian even for complex problems.
So, I tend to believe that your implementation may have some problem. Some 
lessons I learned during my code development:


  *   how do you do the coloring when using PETSc finite differencing? An 
incorrect coloring may give you wrong Jacobian. The simplest way to avoid an 
incorrect coloring is to assume the matrix is dense (slow but error proofing).
  *   Residual function evaluation not correctly implemented can also lead to 
incorrect Jacobian. In your case, you may want to take a careful look at the 
order of execution, when to update your unknown vector, when to perform P1 
reconstruction, and when to evaluate the residual.

-Ling

From: petsc-users  on behalf of Mark Lohry 

Date: Saturday, April 20, 2024 at 1:35 PM
To: PETSc 
Subject: [petsc-users] finite difference jacobian errors when given 
non-constant initial condition
I have a 1-dimensional P1 discontinuous Galerkin discretization of the linear 
advection equation with 4 cells and periodic boundaries on [-pi,+pi]. I'm 
comparing the results from SNESComputeJacobian with a hand-written Jacobian. 
Being linear,
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd
I have a 1-dimensional P1 discontinuous Galerkin discretization of the linear 
advection equation with 4 cells and periodic boundaries on [-pi,+pi]. I'm 
comparing the results from SNESComputeJacobian with a hand-written Jacobian. 
Being linear, the Jacobian should be constant/independent of the solution.

When I set the initial condition passed to SNESComputeJacobian as some 
constant, say f(x)=1 or 0, the petsc finite difference jacobian agrees with my 
hand coded-version. But when I pass it some non-constant value, e.g. 
f(x)=sin(x), something goes horribly wrong in the petsc jacobian. Implementing 
my own rudimentary finite difference approximation (similar to how I thought 
petsc computes it) it returns the correct jacobian to expected error. Any idea 
what could be going on?

Analytically computed Jacobian:
 4.44089e-16 -1.10266  0.31831   -0.085290900   
  -0.31831  1.18795
 1.10266 -4.44089e-16 -1.18795  0.3183100   
 0.0852909 -0.31831
-0.31831  1.18795  4.44089e-16 -1.10266  0.31831   -0.0852909   
 00
   0.0852909 -0.31831  1.10266 -4.44089e-16 -1.18795  0.31831   
 00
   00 -0.31831  1.18795  4.44089e-16 -1.10266   
   0.31831   -0.0852909
   000.0852909 -0.31831  1.10266 -4.44089e-16   
  -1.18795  0.31831
 0.31831   -0.085290900 -0.31831  1.18795  
4.44089e-16 -1.10266
-1.18795  0.31831000.0852909 -0.31831   
   1.10266 -4.44089e-16


petsc finite difference jacobian when given f(x)=1:
 4.44089e-16 -1.10266  0.31831   -0.085290900   
  -0.31831  1.18795
 1.10266 -4.44089e-16 -1.18795  0.3183100   
 0.0852909 -0.31831
-0.31831  1.18795  4.44089e-16 -1.10266  0.31831   -0.0852909   
 00
   0.0852909 -0.31831  1.10266 -4.44089e-16 -1.18795  0.31831   
 00
   00 -0.31831  1.18795  4.44089e-16 -1.10266   
   0.31831   -0.0852909
   000.0852909 -0.31831  1.10266 -4.44089e-16   
  -1.18795  0.31831
 0.31831   -0.085290900 -0.31831  1.18795  
4.44089e-16 -1.10266
-1.18795  0.31831000.0852909 -0.31831   
   1.10266 -4.44089e-16

petsc finite difference jacobian when given f(x) = sin(x):
-1.65547e+08 -3.31856e+08 -1.25427e+09   4.4844e+0800  
1.03206e+08  7.86375e+07
 9.13788e+07  1.83178e+08  6.92336e+08  -2.4753e+0800 
-5.69678e+07 -4.34064e+07
  3.7084e+07  7.43387e+07  2.80969e+08 -1.00455e+08  -5.0384e+07 -2.99747e+07   
 00
  3.7084e+07

Re: [petsc-users] finite difference jacobian errors when given non-constant initial condition

[Zou, Ling via petsc-users]

Hi Mark, I am working on a project having similar numeric you have, 
one-dimensional finite volume method with second-order slope limiter TVD, and 
PETSc finite differencing gives perfect Jacobian even for complex problems.
So, I tend to believe that your implementation may have some problem. Some 
lessons I learned during my code development:


  *   how do you do the coloring when using PETSc finite differencing? An 
incorrect coloring may give you wrong Jacobian. The simplest way to avoid an 
incorrect coloring is to assume the matrix is dense (slow but error proofing).
  *   Residual function evaluation not correctly implemented can also lead to 
incorrect Jacobian. In your case, you may want to take a careful look at the 
order of execution, when to update your unknown vector, when to perform P1 
reconstruction, and when to evaluate the residual.

-Ling

From: petsc-users  on behalf of Mark Lohry 

Date: Saturday, April 20, 2024 at 1:35 PM
To: PETSc 
Subject: [petsc-users] finite difference jacobian errors when given 
non-constant initial condition
I have a 1-dimensional P1 discontinuous Galerkin discretization of the linear 
advection equation with 4 cells and periodic boundaries on [-pi,+pi]. I'm 
comparing the results from SNESComputeJacobian with a hand-written Jacobian. 
Being linear,
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd
I have a 1-dimensional P1 discontinuous Galerkin discretization of the linear 
advection equation with 4 cells and periodic boundaries on [-pi,+pi]. I'm 
comparing the results from SNESComputeJacobian with a hand-written Jacobian. 
Being linear, the Jacobian should be constant/independent of the solution.

When I set the initial condition passed to SNESComputeJacobian as some 
constant, say f(x)=1 or 0, the petsc finite difference jacobian agrees with my 
hand coded-version. But when I pass it some non-constant value, e.g. 
f(x)=sin(x), something goes horribly wrong in the petsc jacobian. Implementing 
my own rudimentary finite difference approximation (similar to how I thought 
petsc computes it) it returns the correct jacobian to expected error. Any idea 
what could be going on?

Analytically computed Jacobian:
 4.44089e-16 -1.10266  0.31831   -0.085290900   
  -0.31831  1.18795
 1.10266 -4.44089e-16 -1.18795  0.3183100   
 0.0852909 -0.31831
-0.31831  1.18795  4.44089e-16 -1.10266  0.31831   -0.0852909   
 00
   0.0852909 -0.31831  1.10266 -4.44089e-16 -1.18795  0.31831   
 00
   00 -0.31831  1.18795  4.44089e-16 -1.10266   
   0.31831   -0.0852909
   000.0852909 -0.31831  1.10266 -4.44089e-16   
  -1.18795  0.31831
 0.31831   -0.085290900 -0.31831  1.18795  
4.44089e-16 -1.10266
-1.18795  0.31831000.0852909 -0.31831   
   1.10266 -4.44089e-16


petsc finite difference jacobian when given f(x)=1:
 4.44089e-16 -1.10266  0.31831   -0.085290900   
  -0.31831  1.18795
 1.10266 -4.44089e-16 -1.18795  0.3183100   
 0.0852909 -0.31831
-0.31831  1.18795  4.44089e-16 -1.10266  0.31831   -0.0852909   
 00
   0.0852909 -0.31831  1.10266 -4.44089e-16 -1.18795  0.31831   
 00
   00 -0.31831  1.18795  4.44089e-16 -1.10266   
   0.31831   -0.0852909
   000.0852909 -0.31831  1.10266 -4.44089e-16   
  -1.18795  0.31831
 0.31831   -0.085290900 -0.31831  1.18795  
4.44089e-16 -1.10266
-1.18795  0.31831000.0852909 -0.31831   
   1.10266 -4.44089e-16

petsc finite difference jacobian when given f(x) = sin(x):
-1.65547e+08 -3.31856e+08 -1.25427e+09   4.4844e+0800  
1.03206e+08  7.86375e+07
 9.13788e+07  1.83178e+08  6.92336e+08  -2.4753e+0800 
-5.69678e+07 -4.34064e+07
  3.7084e+07  7.43387e+07  2.80969e+08 -1.00455e+08  -5.0384e+07 -2.99747e+07   
 00
  3.7084e+07  7.43387e+07  2.80969e+08 -1.00455e+08  -5.0384e+07 -2.99747e+07   
 00
   00  2.80969e+08 -1.00455e+08  -5.0384e+07 -2.99747e+07 
-2.31191e+07 -1.76155e+07
   00  2.80969e+08 -1.00455e+08  -5.0384e+07 -2.99747e+07 
-2.31191e+07 -1.76155e+07
 9.13788e+07  1.83178e+0800 -1.24151e+08 -7.38608e+07 
-5.69678e+07 -4.34064e+07
-1.65547e+08 -3.31856e+0800  2.24919e+08   1.3381e+08  
1.03206e+08  7.86375e+07

Re: [petsc-users] Does ILU(15) still make sense or should just use LU?

[Zou, Ling via petsc-users]

er and you don't 
waste months
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

   You may benefit from a literature search on your model AND preconditioners 
to see what others have used. But I would try PETSc/MUMPS on the biggest size 
you want and see how it goes (better it runs for a little longer and you don't 
waste months trying to find a good preconditioner).





On Mar 28, 2024, at 2:20 PM, Zou, Ling  wrote:

Thank you, Barry.
Yes, I have tried different preconditioners, but in a naïve way, i.e., looping 
through possible options using `-pc_type ` command line.
But no, not in a meaningful way because the lack of understanding of the 
connection between physics (the problem we are dealing with) to math (the 
correct combination of those preconditioners).

-Ling


From: Barry Smith mailto:bsm...@petsc.dev>>
Date: Thursday, March 28, 2024 at 1:09 PM
To: Zou, Ling mailto:l...@anl.gov>>
Cc: Zhang, Hong mailto:hzh...@mcs.anl.gov>>, 
petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov> 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Does ILU(15) still make sense or should just use LU?
1 million is possible for direct solvers using PETSc with the MUMPS direct 
solver when you cannot get a preconditioner to work well for your problems. ILU 
are not very robust preconditioners and I would not rely on them. Have you 
investigated
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

   1 million is possible for direct solvers using PETSc with the MUMPS direct 
solver when you cannot get a preconditioner to work well for your problems.

ILU are not very robust preconditioners and I would not rely on them. Have 
you investigated other preconditioners in PETSc, PCGAMG, PCASM, PCFIELDSPLIT or 
some combination of these preconditioners work for many problems, though 
certainly not all.



On Mar 28, 2024, at 1:14 PM, Zou, Ling mailto:l...@anl.gov>> 
wrote:

Thank you, Barry.
Yeah, this is unfortunate given that the problem we are handling is quite 
heterogeneous (in both mesh and physics).
I expect that our problem sizes will be mostly smaller than 1 million DOF, 
should LU still be a practical solution? Can it scale well if we choose to run 
the problem in a parallel way?

PS1: -ksp_norm_type unpreconditioned did not work as the true residual did not 
go down, even with 300 linear iterations.
PS2: what do you think if it will be beneficial to have more detailed 
discussions (e.g., a presentation?) on the problem we are solving to seek more 
advice?

-Ling


From: Barry Smith mailto:bsm...@petsc.dev>>
Date: Thursday, March 28, 2024 at 11:14 AM
To: Zou, Ling mailto:l...@anl.gov>>
Cc: Zhang, Hong mailto:hzh...@mcs.anl.gov>>, 
petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov> 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Does ILU(15) still make sense or should just use LU?
This is a bad situation, the solver is not really converging. This can happen 
with ILU() sometimes, it so badly scales things that the preconditioned 
residual decreases a lot but the true residual is not really getting smaller. 
Since your matrices
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

   This is a bad situation, the solver is not really converging. This can 
happen with ILU() sometimes, it so badly scales things that the preconditioned 
residual decreases a lot but the true residual is not really getting smaller. 
Since your matrices are small best to stick to LU.

You can use -ksp_norm_type unpreconditioned to force the convergence test 
to use the true residual for a convergence test and the solver will discover 
that it is not converging.

   Barry



On Mar 28, 2024, at 11:43 AM, Zou, Ling via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:

Hong, thanks! That makes perfect sense.
A follow up question about ILU.

The following is the performance of ILU(5). Note that each KPS solving reports 
converged but as the output shows, the preconditioned residual does while true 
residual does not. Is there any way this performance could be improved?
Background: the preconditioning matrix is finite difference generated, and 
should be exact.

-Ling

Time Step 21, time = -491.75, dt = 1
NL Step =  0, fnorm =  6.98749E+01
0 KSP preconditioned resid norm 1.684131526824e+04 true resid norm 
6.987489798042e+01 ||r(i)||/||b|| 1.e+00
1 KSP preconditioned resid norm 5.970568556551e+02 true resid norm 
6.459553545222e+01 ||r(i)||/||b|| 9.244455064582e-01
2 KSP preconditioned resid norm 3.349113985192e+02 true resid norm 
7.250836872274e+01 ||r(i)||/||b|| 1.037688366186e+00
3 KSP preconditioned resid norm 3.290585904777e+01 true resid norm 
1.186282435163e+

Re: [petsc-users] Does ILU(15) still make sense or should just use LU?

[Zou, Ling via petsc-users]

rely on them. Have you 
investigated

ZjQcmQRYFpfptBannerStart

This Message Is From an External Sender

This message came from outside your organization.



ZjQcmQRYFpfptBannerEnd



   1 million is possible for direct solvers using PETSc with the MUMPS direct 
solver when you cannot get a preconditioner to work well for your problems.



ILU are not very robust preconditioners and I would not rely on them. Have 
you investigated other preconditioners in PETSc, PCGAMG, PCASM, PCFIELDSPLIT or 
some combination of these preconditioners work for many problems, though 
certainly not all.





On Mar 28, 2024, at 1:14 PM, Zou, Ling mailto:l...@anl.gov>> 
wrote:



Thank you, Barry.

Yeah, this is unfortunate given that the problem we are handling is quite 
heterogeneous (in both mesh and physics).

I expect that our problem sizes will be mostly smaller than 1 million DOF, 
should LU still be a practical solution? Can it scale well if we choose to run 
the problem in a parallel way?



PS1: -ksp_norm_type unpreconditioned did not work as the true residual did not 
go down, even with 300 linear iterations.

PS2: what do you think if it will be beneficial to have more detailed 
discussions (e.g., a presentation?) on the problem we are solving to seek more 
advice?



-Ling



From: Barry Smith mailto:bsm...@petsc.dev>>
Date: Thursday, March 28, 2024 at 11:14 AM
To: Zou, Ling mailto:l...@anl.gov>>
Cc: Zhang, Hong mailto:hzh...@mcs.anl.gov>>, 
petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov> 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Does ILU(15) still make sense or should just use LU?

This is a bad situation, the solver is not really converging. This can happen 
with ILU() sometimes, it so badly scales things that the preconditioned 
residual decreases a lot but the true residual is not really getting smaller. 
Since your matrices

ZjQcmQRYFpfptBannerStart

This Message Is From an External Sender

This message came from outside your organization.



ZjQcmQRYFpfptBannerEnd



   This is a bad situation, the solver is not really converging. This can 
happen with ILU() sometimes, it so badly scales things that the preconditioned 
residual decreases a lot but the true residual is not really getting smaller. 
Since your matrices are small best to stick to LU.



You can use -ksp_norm_type unpreconditioned to force the convergence test 
to use the true residual for a convergence test and the solver will discover 
that it is not converging.



   Barry





On Mar 28, 2024, at 11:43 AM, Zou, Ling via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:



Hong, thanks! That makes perfect sense.

A follow up question about ILU.



The following is the performance of ILU(5). Note that each KPS solving reports 
converged but as the output shows, the preconditioned residual does while true 
residual does not. Is there any way this performance could be improved?

Background: the preconditioning matrix is finite difference generated, and 
should be exact.



-Ling



Time Step 21, time = -491.75, dt = 1

NL Step =  0, fnorm =  6.98749E+01

0 KSP preconditioned resid norm 1.684131526824e+04 true resid norm 
6.987489798042e+01 ||r(i)||/||b|| 1.e+00

1 KSP preconditioned resid norm 5.970568556551e+02 true resid norm 
6.459553545222e+01 ||r(i)||/||b|| 9.244455064582e-01

2 KSP preconditioned resid norm 3.349113985192e+02 true resid norm 
7.250836872274e+01 ||r(i)||/||b|| 1.037688366186e+00

3 KSP preconditioned resid norm 3.290585904777e+01 true resid norm 
1.186282435163e+02 ||r(i)||/||b|| 1.697723316169e+00

4 KSP preconditioned resid norm 8.530606201233e+00 true resid norm 
4.088729421459e+01 ||r(i)||/||b|| 5.851499665310e-01

  Linear solve converged due to CONVERGED_RTOL iterations 4

NL Step =  1, fnorm =  4.08788E+01

0 KSP preconditioned resid norm 1.851047973094e+03 true resid norm 
4.087882723223e+01 ||r(i)||/||b|| 1.e+00

1 KSP preconditioned resid norm 3.696809614513e+01 true resid norm 
2.720016413105e+01 ||r(i)||/||b|| 6.653851387793e-01

2 KSP preconditioned resid norm 5.751891392534e+00 true resid norm 
3.326338240872e+01 ||r(i)||/||b|| 8.137068663873e-01

3 KSP preconditioned resid norm 8.540729397958e-01 true resid norm 
8.672410748720e+00 ||r(i)||/||b|| 2.121492062249e-01

  Linear solve converged due to CONVERGED_RTOL iterations 3

NL Step =  2, fnorm =  8.67124E+00

0 KSP preconditioned resid norm 5.511333966852e+00 true resid norm 
8.671237519593e+00 ||r(i)||/||b|| 1.e+00

1 KSP preconditioned resid norm 1.174962622023e+00 true resid norm 
8.731034658309e+00 ||r(i)||/||b|| 1.006896032842e+00

2 KSP preconditioned resid norm 1.104604471016e+00 true resid norm 
1.018397505468e+01 ||r(i)||/||b|| 1.174454630227e+00

3 KSP preconditioned resid norm 4.257063674222e-01 true resid norm 
4.023093124996e+00 ||r(i)||/||b|| 4.639583584126e-01

4 KSP

Re: [petsc-users] Does ILU(15) still make sense or should just use LU?

[Zou, Ling via petsc-users]

Thank you. Those are great suggestions. Although I mentioned 1 million DOF, but 
we rarely go there, so maybe stick with what is working now, and meanwhile 
seeking helps from literatures.

-Ling

From: Barry Smith 
Date: Thursday, March 28, 2024 at 2:26 PM
To: Zou, Ling 
Cc: Zhang, Hong , petsc-users@mcs.anl.gov 

Subject: Re: [petsc-users] Does ILU(15) still make sense or should just use LU?
You may benefit from a literature search on your model AND preconditioners to 
see what others have used. But I would try PETSc/MUMPS on the biggest size you 
want and see how it goes (better it runs for a little longer and you don't 
waste months
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

   You may benefit from a literature search on your model AND preconditioners 
to see what others have used. But I would try PETSc/MUMPS on the biggest size 
you want and see how it goes (better it runs for a little longer and you don't 
waste months trying to find a good preconditioner).





On Mar 28, 2024, at 2:20 PM, Zou, Ling  wrote:

Thank you, Barry.
Yes, I have tried different preconditioners, but in a naïve way, i.e., looping 
through possible options using `-pc_type ` command line.
But no, not in a meaningful way because the lack of understanding of the 
connection between physics (the problem we are dealing with) to math (the 
correct combination of those preconditioners).

-Ling

From: Barry Smith mailto:bsm...@petsc.dev>>
Date: Thursday, March 28, 2024 at 1:09 PM
To: Zou, Ling mailto:l...@anl.gov>>
Cc: Zhang, Hong mailto:hzh...@mcs.anl.gov>>, 
petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov> 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Does ILU(15) still make sense or should just use LU?
1 million is possible for direct solvers using PETSc with the MUMPS direct 
solver when you cannot get a preconditioner to work well for your problems. ILU 
are not very robust preconditioners and I would not rely on them. Have you 
investigated
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

   1 million is possible for direct solvers using PETSc with the MUMPS direct 
solver when you cannot get a preconditioner to work well for your problems.

ILU are not very robust preconditioners and I would not rely on them. Have 
you investigated other preconditioners in PETSc, PCGAMG, PCASM, PCFIELDSPLIT or 
some combination of these preconditioners work for many problems, though 
certainly not all.


On Mar 28, 2024, at 1:14 PM, Zou, Ling mailto:l...@anl.gov>> 
wrote:

Thank you, Barry.
Yeah, this is unfortunate given that the problem we are handling is quite 
heterogeneous (in both mesh and physics).
I expect that our problem sizes will be mostly smaller than 1 million DOF, 
should LU still be a practical solution? Can it scale well if we choose to run 
the problem in a parallel way?

PS1: -ksp_norm_type unpreconditioned did not work as the true residual did not 
go down, even with 300 linear iterations.
PS2: what do you think if it will be beneficial to have more detailed 
discussions (e.g., a presentation?) on the problem we are solving to seek more 
advice?

-Ling

From: Barry Smith mailto:bsm...@petsc.dev>>
Date: Thursday, March 28, 2024 at 11:14 AM
To: Zou, Ling mailto:l...@anl.gov>>
Cc: Zhang, Hong mailto:hzh...@mcs.anl.gov>>, 
petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov> 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Does ILU(15) still make sense or should just use LU?
This is a bad situation, the solver is not really converging. This can happen 
with ILU() sometimes, it so badly scales things that the preconditioned 
residual decreases a lot but the true residual is not really getting smaller. 
Since your matrices
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

   This is a bad situation, the solver is not really converging. This can 
happen with ILU() sometimes, it so badly scales things that the preconditioned 
residual decreases a lot but the true residual is not really getting smaller. 
Since your matrices are small best to stick to LU.

You can use -ksp_norm_type unpreconditioned to force the convergence test 
to use the true residual for a convergence test and the solver will discover 
that it is not converging.

   Barry


On Mar 28, 2024, at 11:43 AM, Zou, Ling via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:

Hong, thanks! That makes perfect sense.
A follow up question about ILU.

The following is the performance of ILU(5). Note that each KPS solving reports 
converged but as the output shows, the preconditioned residual does while true 
residual does not. Is there any way this performance could be improve

Re: [petsc-users] Does ILU(15) still make sense or should just use LU?

[Zou, Ling via petsc-users]

Thank you, Barry.
Yes, I have tried different preconditioners, but in a naïve way, i.e., looping 
through possible options using `-pc_type ` command line.
But no, not in a meaningful way because the lack of understanding of the 
connection between physics (the problem we are dealing with) to math (the 
correct combination of those preconditioners).

-Ling

From: Barry Smith 
Date: Thursday, March 28, 2024 at 1:09 PM
To: Zou, Ling 
Cc: Zhang, Hong , petsc-users@mcs.anl.gov 

Subject: Re: [petsc-users] Does ILU(15) still make sense or should just use LU?
1 million is possible for direct solvers using PETSc with the MUMPS direct 
solver when you cannot get a preconditioner to work well for your problems. ILU 
are not very robust preconditioners and I would not rely on them. Have you 
investigated
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

   1 million is possible for direct solvers using PETSc with the MUMPS direct 
solver when you cannot get a preconditioner to work well for your problems.

ILU are not very robust preconditioners and I would not rely on them. Have 
you investigated other preconditioners in PETSc, PCGAMG, PCASM, PCFIELDSPLIT or 
some combination of these preconditioners work for many problems, though 
certainly not all.



On Mar 28, 2024, at 1:14 PM, Zou, Ling  wrote:

Thank you, Barry.
Yeah, this is unfortunate given that the problem we are handling is quite 
heterogeneous (in both mesh and physics).
I expect that our problem sizes will be mostly smaller than 1 million DOF, 
should LU still be a practical solution? Can it scale well if we choose to run 
the problem in a parallel way?

PS1: -ksp_norm_type unpreconditioned did not work as the true residual did not 
go down, even with 300 linear iterations.
PS2: what do you think if it will be beneficial to have more detailed 
discussions (e.g., a presentation?) on the problem we are solving to seek more 
advice?

-Ling

From: Barry Smith mailto:bsm...@petsc.dev>>
Date: Thursday, March 28, 2024 at 11:14 AM
To: Zou, Ling mailto:l...@anl.gov>>
Cc: Zhang, Hong mailto:hzh...@mcs.anl.gov>>, 
petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov> 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Does ILU(15) still make sense or should just use LU?
This is a bad situation, the solver is not really converging. This can happen 
with ILU() sometimes, it so badly scales things that the preconditioned 
residual decreases a lot but the true residual is not really getting smaller. 
Since your matrices
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

   This is a bad situation, the solver is not really converging. This can 
happen with ILU() sometimes, it so badly scales things that the preconditioned 
residual decreases a lot but the true residual is not really getting smaller. 
Since your matrices are small best to stick to LU.

You can use -ksp_norm_type unpreconditioned to force the convergence test 
to use the true residual for a convergence test and the solver will discover 
that it is not converging.

   Barry


On Mar 28, 2024, at 11:43 AM, Zou, Ling via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:

Hong, thanks! That makes perfect sense.
A follow up question about ILU.

The following is the performance of ILU(5). Note that each KPS solving reports 
converged but as the output shows, the preconditioned residual does while true 
residual does not. Is there any way this performance could be improved?
Background: the preconditioning matrix is finite difference generated, and 
should be exact.

-Ling

Time Step 21, time = -491.75, dt = 1
NL Step =  0, fnorm =  6.98749E+01
0 KSP preconditioned resid norm 1.684131526824e+04 true resid norm 
6.987489798042e+01 ||r(i)||/||b|| 1.e+00
1 KSP preconditioned resid norm 5.970568556551e+02 true resid norm 
6.459553545222e+01 ||r(i)||/||b|| 9.244455064582e-01
2 KSP preconditioned resid norm 3.349113985192e+02 true resid norm 
7.250836872274e+01 ||r(i)||/||b|| 1.037688366186e+00
3 KSP preconditioned resid norm 3.290585904777e+01 true resid norm 
1.186282435163e+02 ||r(i)||/||b|| 1.697723316169e+00
4 KSP preconditioned resid norm 8.530606201233e+00 true resid norm 
4.088729421459e+01 ||r(i)||/||b|| 5.851499665310e-01
  Linear solve converged due to CONVERGED_RTOL iterations 4
NL Step =  1, fnorm =  4.08788E+01
0 KSP preconditioned resid norm 1.851047973094e+03 true resid norm 
4.087882723223e+01 ||r(i)||/||b|| 1.e+00
1 KSP preconditioned resid norm 3.696809614513e+01 true resid norm 
2.720016413105e+01 ||r(i)||/||b|| 6.653851387793e-01
2 KSP preconditioned resid norm 5.751891392534e+00 true resid norm 
3.326338240872e+01 ||r(i)||/||b|| 8.137068663873e-01
3 KSP preconditioned resid norm 8.540729397958e

Re: [petsc-users] Does ILU(15) still make sense or should just use LU?

[Zou, Ling via petsc-users]

Thank you, Barry.
Yeah, this is unfortunate given that the problem we are handling is quite 
heterogeneous (in both mesh and physics).
I expect that our problem sizes will be mostly smaller than 1 million DOF, 
should LU still be a practical solution? Can it scale well if we choose to run 
the problem in a parallel way?

PS1: -ksp_norm_type unpreconditioned did not work as the true residual did not 
go down, even with 300 linear iterations.
PS2: what do you think if it will be beneficial to have more detailed 
discussions (e.g., a presentation?) on the problem we are solving to seek more 
advice?

-Ling

From: Barry Smith 
Date: Thursday, March 28, 2024 at 11:14 AM
To: Zou, Ling 
Cc: Zhang, Hong , petsc-users@mcs.anl.gov 

Subject: Re: [petsc-users] Does ILU(15) still make sense or should just use LU?
This is a bad situation, the solver is not really converging. This can happen 
with ILU() sometimes, it so badly scales things that the preconditioned 
residual decreases a lot but the true residual is not really getting smaller. 
Since your matrices
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

   This is a bad situation, the solver is not really converging. This can 
happen with ILU() sometimes, it so badly scales things that the preconditioned 
residual decreases a lot but the true residual is not really getting smaller. 
Since your matrices are small best to stick to LU.

You can use -ksp_norm_type unpreconditioned to force the convergence test 
to use the true residual for a convergence test and the solver will discover 
that it is not converging.

   Barry



On Mar 28, 2024, at 11:43 AM, Zou, Ling via petsc-users 
 wrote:

Hong, thanks! That makes perfect sense.
A follow up question about ILU.

The following is the performance of ILU(5). Note that each KPS solving reports 
converged but as the output shows, the preconditioned residual does while true 
residual does not. Is there any way this performance could be improved?
Background: the preconditioning matrix is finite difference generated, and 
should be exact.

-Ling

Time Step 21, time = -491.75, dt = 1
NL Step =  0, fnorm =  6.98749E+01
0 KSP preconditioned resid norm 1.684131526824e+04 true resid norm 
6.987489798042e+01 ||r(i)||/||b|| 1.e+00
1 KSP preconditioned resid norm 5.970568556551e+02 true resid norm 
6.459553545222e+01 ||r(i)||/||b|| 9.244455064582e-01
2 KSP preconditioned resid norm 3.349113985192e+02 true resid norm 
7.250836872274e+01 ||r(i)||/||b|| 1.037688366186e+00
3 KSP preconditioned resid norm 3.290585904777e+01 true resid norm 
1.186282435163e+02 ||r(i)||/||b|| 1.697723316169e+00
4 KSP preconditioned resid norm 8.530606201233e+00 true resid norm 
4.088729421459e+01 ||r(i)||/||b|| 5.851499665310e-01
  Linear solve converged due to CONVERGED_RTOL iterations 4
NL Step =  1, fnorm =  4.08788E+01
0 KSP preconditioned resid norm 1.851047973094e+03 true resid norm 
4.087882723223e+01 ||r(i)||/||b|| 1.e+00
1 KSP preconditioned resid norm 3.696809614513e+01 true resid norm 
2.720016413105e+01 ||r(i)||/||b|| 6.653851387793e-01
2 KSP preconditioned resid norm 5.751891392534e+00 true resid norm 
3.326338240872e+01 ||r(i)||/||b|| 8.137068663873e-01
3 KSP preconditioned resid norm 8.540729397958e-01 true resid norm 
8.672410748720e+00 ||r(i)||/||b|| 2.121492062249e-01
  Linear solve converged due to CONVERGED_RTOL iterations 3
NL Step =  2, fnorm =  8.67124E+00
0 KSP preconditioned resid norm 5.511333966852e+00 true resid norm 
8.671237519593e+00 ||r(i)||/||b|| 1.e+00
1 KSP preconditioned resid norm 1.174962622023e+00 true resid norm 
8.731034658309e+00 ||r(i)||/||b|| 1.006896032842e+00
2 KSP preconditioned resid norm 1.104604471016e+00 true resid norm 
1.018397505468e+01 ||r(i)||/||b|| 1.174454630227e+00
3 KSP preconditioned resid norm 4.257063674222e-01 true resid norm 
4.023093124996e+00 ||r(i)||/||b|| 4.639583584126e-01
4 KSP preconditioned resid norm 1.023038868263e-01 true resid norm 
2.365298462869e+00 ||r(i)||/||b|| 2.727751901068e-01
5 KSP preconditioned resid norm 4.073772638935e-02 true resid norm 
2.302623112025e+00 ||r(i)||/||b|| 2.655472309255e-01
6 KSP preconditioned resid norm 1.510323179379e-02 true resid norm 
2.300216593521e+00 ||r(i)||/||b|| 2.652697020839e-01
7 KSP preconditioned resid norm 1.337324816903e-02 true resid norm 
2.300057733345e+00 ||r(i)||/||b|| 2.652513817259e-01
8 KSP preconditioned resid norm 1.247384902656e-02 true resid norm 
2.300456226062e+00 ||r(i)||/||b|| 2.652973374174e-01
9 KSP preconditioned resid norm 1.247038855375e-02 true resid norm 
2.300532560993e+00 ||r(i)||/||b|| 2.653061406512e-01
   10 KSP preconditioned resid norm 1.244611343317e-02 true resid norm 
2.299441241514e+00 ||r(i)||/||b|| 2.651802855496e-01
   11 KSP preconditioned resid norm 1.227243209527e-02 true resid

Re: [petsc-users] Does ILU(15) still make sense or should just use LU?

[Zou, Ling via petsc-users]

, March 27, 2024 at 4:59 PM
To: petsc-users@mcs.anl.gov , Zou, Ling 
Subject: Re: Does ILU(15) still make sense or should just use LU?
Ling,
ILU(level) is used for saving storage space with more computations. Normally, 
we use level=1 or 2. It does not make sense to use level 15. If you have 
sufficient space, LU would be the best.
Hong


From: petsc-users  on behalf of Zou, Ling via 
petsc-users 
Sent: Wednesday, March 27, 2024 4:24 PM
To: petsc-users@mcs.anl.gov 
Subject: [petsc-users] Does ILU(15) still make sense or should just use LU?


Hi, I’d like to avoid using LU, but in some cases to use ILU and still 
converge, I have to go to ILU(15), i.e., `-pc_factor_levels 15`. Does it still 
make sense, or should I give it up and switch to LU?



For this particular case, ~2k DoF, and both ILU(15) and LU perform similarly in 
terms of wall time.



-Ling

[petsc-users] Does ILU(15) still make sense or should just use LU?

[Zou, Ling via petsc-users]

Hi, I’d like to avoid using LU, but in some cases to use ILU and still 
converge, I have to go to ILU(15), i.e., `-pc_factor_levels 15`. Does it still 
make sense, or should I give it up and switch to LU?

For this particular case, ~2k DoF, and both ILU(15) and LU perform similarly in 
terms of wall time.

-Ling

[petsc-users] Possible bug associated with '-pc_factor_nonzeros_along_diagonal' option?

[Zou, Ling via petsc-users]

Dear all,

For the same code and same input file, using 
'-pc_factor_nonzeros_along_diagonal' causing a seg fault, while 
'-pc_factor_shift_type nonzero -pc_factor_shift_amount 1e-8' works fine.
The seg fault happened after the very first residual vec norm was evaluated, 
somewhere before stepping into the second residual vec norm evaluation, i.e.,


NL Step =  0, fnorm =  5.22084E+01

Segmentation fault: 11

I tend to believe that it is not in my `snesFormFunction` because I used print 
to make sure to see that the ‘sef fault’ is after snesFormFunction was called 
in the first residual vec norm was shown, and before snesFormFunction was 
called for the next residual vec norm.

I am behind ‘moose,’ so the debug version did not tell much info, showing the 
exact error message.
It’s not an urgent issue for me, since the other option worked fine.

Best,

-Ling

Re: [petsc-users] 'Preconditioning' with lower-order method

[Zou, Ling via petsc-users]

From: Zhang, Hong 
Date: Monday, March 4, 2024 at 6:34 PM
To: Zou, Ling 
Cc: Jed Brown , Barry Smith , 
petsc-users@mcs.anl.gov 
Subject: Re: [petsc-users] 'Preconditioning' with lower-order method
Ling,

Are you using PETSc TS? If so, it may worth trying Crank-Nicolson first to see 
if the nonlinear solve becomes faster.

>>>>> No, I’m not using TS, and I don’t plan to use CN. From my experience, 
>>>>> when dealing with (nearly) incompressible flow problems, CN often cause 
>>>>> (very large) pressure temporal oscillations, and to avoid that, the 
>>>>> pressure is often using fully implicit method, so that would cause quite 
>>>>> some code implementation issue.
For the pressure oscillation issue, also see page 7 of INL/EXT-12-27197. Notes 
on Newton-Krylov Based Incompressible Flow Projection Solver


In addition, you can try to improve the performance by pruning the Jacobian 
matrix.

TSPruneIJacobianColor()<https://urldefense.us/v3/__https://petsc.org/main/manualpages/TS/TSPruneIJacobianColor/__;!!G_uCfscf7eWS!aXVa0uz1LUIOdvZEPlRJOhRzz9h8MSM4vhl93kknxKGb8hkTyjCFJmSZIGr0fYx90rrqotBGdw-NwiF27Ec$
 > sometimes can reduce the number of colors especially for high-order methods 
and make your Jacobian matrix more compact. An example of usage can be found 
here<https://urldefense.us/v3/__https://petsc.org/release/src/ts/tests/ex5.c.html__;!!G_uCfscf7eWS!aXVa0uz1LUIOdvZEPlRJOhRzz9h8MSM4vhl93kknxKGb8hkTyjCFJmSZIGr0fYx90rrqotBGdw-NmisGfPQ$
 >. If you are not using TS, there is a SNES version 
SNESPruneJacobianColor()<https://urldefense.us/v3/__https://petsc.org/main/manualpages/SNES/SNESPruneJacobianColor/__;!!G_uCfscf7eWS!aXVa0uz1LUIOdvZEPlRJOhRzz9h8MSM4vhl93kknxKGb8hkTyjCFJmSZIGr0fYx90rrqotBGdw-NfYBUJqI$
 > for the same functionality.

>>>>> The following code is how I setup the coloring.
  {
// Create Matrix-free context
MatCreateSNESMF(snes, _MatrixFree);

// Let the problem setup Jacobian matrix sparsity
p_sim->FillJacobianMatrixNonZeroEntry(P_Mat);

// See PETSc examples:
// 
https://urldefense.us/v3/__https://petsc.org/release/src/snes/tutorials/ex14.c.html__;!!G_uCfscf7eWS!aXVa0uz1LUIOdvZEPlRJOhRzz9h8MSM4vhl93kknxKGb8hkTyjCFJmSZIGr0fYx90rrqotBGdw-N3ZHE6Qw$
 
// 
https://urldefense.us/v3/__https://petsc.org/release/src/mat/tutorials/ex16.c.html__;!!G_uCfscf7eWS!aXVa0uz1LUIOdvZEPlRJOhRzz9h8MSM4vhl93kknxKGb8hkTyjCFJmSZIGr0fYx90rrqotBGdw-NcpSrhMM$
 
ISColoring iscoloring;
MatColoring mc;
MatColoringCreate(P_Mat, );
MatColoringSetType(mc, MATCOLORINGSL);
MatColoringSetFromOptions(mc);
MatColoringApply(mc, );
MatColoringDestroy();
MatFDColoringCreate(P_Mat, iscoloring, );
MatFDColoringSetFunction(
fdcoloring, (PetscErrorCode(*)(void))(void (*)(void))snesFormFunction, 
this);
MatFDColoringSetFromOptions(fdcoloring);
MatFDColoringSetUp(P_Mat, iscoloring, fdcoloring);
ISColoringDestroy();

// Should I prune here? Like

SNESPruneJacobianColor<https://urldefense.us/v3/__https://petsc.org/main/manualpages/SNES/SNESPruneJacobianColor/__;!!G_uCfscf7eWS!aXVa0uz1LUIOdvZEPlRJOhRzz9h8MSM4vhl93kknxKGb8hkTyjCFJmSZIGr0fYx90rrqotBGdw-NfYBUJqI$
 >(snes, P_Mat, P_Mat);


SNESSetJacobian(snes,// snes
J_MatrixFree,// Jacobian-free
P_Mat,   // Preconditioning matrix
SNESComputeJacobianDefaultColor, // Use finite differencing 
and coloring
    fdcoloring);     // fdcoloring
  }

Thanks,

-Ling




Hong (Mr.)


On Mar 3, 2024, at 11:48 PM, Zou, Ling via petsc-users 
 wrote:



From: Jed Brown mailto:j...@jedbrown.org>>
Date: Sunday, March 3, 2024 at 11:35 PM
To: Zou, Ling mailto:l...@anl.gov>>, Barry Smith 
mailto:bsm...@petsc.dev>>
Cc: petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov> 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] 'Preconditioning' with lower-order method
If you're having PETSc use coloring and have confirmed that the stencil is 
sufficient, then it would be nonsmoothness (again, consider the limiter you've 
chosen) preventing quadratic convergence (assuming that doesn't kick in 
eventually). Note
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

If you're having PETSc use coloring and have confirmed that the stencil is 
sufficient, then it would be nonsmoothness (again, consider the limiter you've 
chosen) preventing quadratic convergence (assuming that doesn't kick in 
eventually).



• Yes, I do use coloring, and I do provide sufficient stencil, i.e., 
neighbor’s neighbor. The sufficiency is confirmed by PETSc’s 
-snes_test_jacobian and -snes_test_jacobian_view options.



Note that assem

Re: [petsc-users] 'Preconditioning' with lower-order method

[Zou, Ling via petsc-users]

 Object: 1 MPI process

  type: newtonls

  maximum iterations=20, maximum function evaluations=1

  tolerances: relative=1e-08, absolute=1e-06, solution=1e-08

  total number of linear solver iterations=44

  total number of function evaluations=51

  norm schedule ALWAYS

  Jacobian is applied matrix-free with differencing

  Preconditioning Jacobian is built using finite differences with coloring

  SNESLineSearch Object: 1 MPI process

type: bt

  interpolation: cubic

  alpha=1.00e-04

maxstep=1.00e+08, minlambda=1.00e-12

tolerances: relative=1.00e-08, absolute=1.00e-15, 
lambda=1.00e-08

maximum iterations=40

  KSP Object: 1 MPI process

type: gmres

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

  happy breakdown tolerance 1e-30

maximum iterations=100, initial guess is zero

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

left preconditioning

using PRECONDITIONED norm type for convergence test

  PC Object: 1 MPI process

type: ilu

  out-of-place factorization

  0 levels of fill

  tolerance for zero pivot 2.22045e-14

  using diagonal shift to prevent zero pivot [NONZERO]

  matrix ordering: rcm

  factor fill ratio given 1., needed 1.

Factored matrix follows:

  Mat Object: 1 MPI process

type: seqaij

rows=8715, cols=8715

package used to perform factorization: petsc

total: nonzeros=38485, allocated nonzeros=38485

  not using I-node routines

linear system matrix followed by preconditioner matrix:

Mat Object: 1 MPI process

  type: mffd

  rows=8715, cols=8715

Matrix-free approximation:

  err=1.49012e-08 (relative error in function evaluation)

  Using wp compute h routine

  Does not compute normU

Mat Object: 1 MPI process

  type: seqaij

  rows=8715, cols=8715

  total: nonzeros=38485, allocated nonzeros=38485

  total number of mallocs used during MatSetValues calls=0

not using I-node routines

 Solve Converged!







"Zou, Ling via petsc-users"  writes:



> Barry, thank you.

> I am not sure if I exactly follow you on this:

> “Are you forming the Jacobian for the first and second order cases inside of 
> Newton?”

>

> The problem that we deal with, heat/mass transfer in heterogeneous systems 
> (reactor system), is generally small in terms of size, i.e., # of DOFs 
> (several k to maybe 100k level), so for now, I completely rely on PETSc to 
> compute Jacobian, i.e., finite-differencing.

>

> That’s a good suggestion to see the time spent during various events.

> What motivated me to try the options are the following observations.

>

> 2nd order FVM:

>

> Time Step 149, time = 13229.7, dt = 100

>

> NL Step =  0, fnorm =  7.80968E-03

>

> NL Step =  1, fnorm =  7.65731E-03

>

> NL Step =  2, fnorm =  6.85034E-03

>

> NL Step =  3, fnorm =  6.11873E-03

>

> NL Step =  4, fnorm =  1.57347E-03

>

> NL Step =  5, fnorm =  9.03536E-04

>

>  Solve Converged!

>

> 1st order FVM:

>

> Time Step 149, time = 13229.7, dt = 100

>

> NL Step =  0, fnorm =  7.90072E-03

>

> NL Step =  1, fnorm =  2.01919E-04

>

> NL Step =  2, fnorm =  1.06960E-05

>

> NL Step =  3, fnorm =  2.41683E-09

>

>  Solve Converged!

>

> Notice the obvious ‘stagnant’ in residual for the 2nd order method while not 
> in the 1st order.

> For the same problem, the wall time is 10 sec vs 6 sec. I would be happy if I 
> can reduce 2 sec for the 2nd order method.

>

> -Ling

>

> From: Barry Smith 

> Date: Sunday, March 3, 2024 at 12:06 PM

> To: Zou, Ling 

> Cc: petsc-users@mcs.anl.gov 

> Subject: Re: [petsc-users] 'Preconditioning' with lower-order method

> Are you forming the Jacobian for the first and second order cases inside of 
> Newton? You can run both with -log_view to see how much time is spent in the 
> various events (compute function, compute Jacobian, linear solve, .. . ) for 
> the two cases

> ZjQcmQRYFpfptBannerStart

> This Message Is From an External Sender

> This message came from outside your organization.

>

> ZjQcmQRYFpfptBannerEnd

>

>Are you forming the Jacobian for the first and second order cases inside 
> of Newton?

>

>You can run both with -log_view to see how much time is spent in the 
> various events (compute function, compute Jacobian, linear solve, ...) for 
> the two cases and compare them.

>

>

>

>

> On Mar 3, 2024, at 11:42 AM, Zou, Ling via petsc-users 
>  wrote:

>

> Original email may have been sent to the inc

Re: [petsc-users] FW: 'Preconditioning' with lower-order method

[Zou, Ling via petsc-users]

From: Jed Brown 
Date: Sunday, March 3, 2024 at 9:24 PM
To: Zou, Ling , petsc-users@mcs.anl.gov 
Subject: Re: [petsc-users] FW: 'Preconditioning' with lower-order method
One option is to form the preconditioner using the FV1 method, which is sparser 
and satisfies h-ellipticity, while using FV2 for the residual and (optionally) 
for matrix-free operator application. FV1 is a highly diffusive method so in a 
sense,
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

One option is to form the preconditioner using the FV1 method, which is sparser 
and satisfies h-ellipticity, while using FV2 for the residual and (optionally) 
for matrix-free operator application.



<<< In terms of code implementation, this seems a bit tricky to me. Looks to me 
that I have to know exactly who is calling the residual function, if its MF 
operation, using FV2, while if finite-differencing for Jacobian, using FV1. 
Currently, I don’t know how to do it.

Another thing I’d like to mention is that the linear solver has never really 
been an issue. While the non-linear solver for the FV2 scheme often ‘stagnant’ 
at the first a couple of non-linear iteration [see the other email reply to 
Barry]. Seems to me, the additional nonlinearity from the TVD limiter causing 
difficulty to PETSc to find the attraction zone.



FV1 is a highly diffusive method so in a sense, it's much less faithful to the 
physics and (say, in the case of fluids) similar to a much lower-Reynolds 
number (if you use a modified equation analysis to work out the effective 
Reynolds number in the presence of the numerical diffusion).



It's good to put some thought into your choice of limiter. Note that 
intersection of second order and TVD methods leads to mandatory nonsmoothness 
(discontinuous derivatives).



<<< Yeah… I am afraid that the TVD limiter is the issue, so that’s the reason 
I’d try to use FV1 to bring the solution (hopefully) closer to the real 
solution so the non-linear solver has an easy job to do.



"Zou, Ling via petsc-users"  writes:



> Original email may have been sent to the incorrect place.

> See below.

>

> -Ling

>

> From: Zou, Ling 

> Date: Sunday, March 3, 2024 at 10:34 AM

> To: petsc-users 

> Subject: 'Preconditioning' with lower-order method

> Hi all,

>

> I am solving a PDE system over a spatial domain. Numerical methods are:

>

>   *   Finite volume method (both 1st and 2nd order implemented)

>   *   BDF1 and BDF2 for time integration.

> What I have noticed is that 1st order FVM converges much faster than 2nd 
> order FVM, regardless the time integration scheme. Well, not surprising since 
> 2nd order FVM introduces additional non-linearity.

>

> I’m thinking about two possible ways to speed up 2nd order FVM, and would 
> like to get some thoughts or community knowledge before jumping into code 
> implementation.

>

> Say, let the 2nd order FVM residual function be F2(x) = 0; and the 1st order 
> FVM residual function be F1(x) = 0.

>

>   1.  Option – 1, multi-step for each time step

> Step 1: solving F1(x) = 0 to obtain a temporary solution x1

> Step 2: feed x1 as an initial guess to solve F2(x) = 0 to obtain the final 
> solution.

> [Not sure if gain any saving at all]

>

>

>   1.  Option -2, dynamically changing residual function F(x)

>

> In pseudo code, would be something like.

>

>

>

> snesFormFunction(SNES snes, Vec u, Vec f, void *)

>

> {

>

>   if (snes.nl_it_no < 4) // 4 being arbitrary here

>

> f = F1(u);

>

>   else

>

> f = F2(u);

>

> }

>

>

>

> I know this might be a bit crazy since it may crash after switching residual 
> function, still, any thoughts?

>

> Best,

>

> -Ling

Re: [petsc-users] 'Preconditioning' with lower-order method

[Zou, Ling via petsc-users]

Barry, thank you.
I am not sure if I exactly follow you on this:
“Are you forming the Jacobian for the first and second order cases inside of 
Newton?”

The problem that we deal with, heat/mass transfer in heterogeneous systems 
(reactor system), is generally small in terms of size, i.e., # of DOFs (several 
k to maybe 100k level), so for now, I completely rely on PETSc to compute 
Jacobian, i.e., finite-differencing.

That’s a good suggestion to see the time spent during various events.
What motivated me to try the options are the following observations.

2nd order FVM:

Time Step 149, time = 13229.7, dt = 100

NL Step =  0, fnorm =  7.80968E-03

NL Step =  1, fnorm =  7.65731E-03

NL Step =  2, fnorm =  6.85034E-03

NL Step =  3, fnorm =  6.11873E-03

NL Step =  4, fnorm =  1.57347E-03

NL Step =  5, fnorm =  9.03536E-04

 Solve Converged!

1st order FVM:

Time Step 149, time = 13229.7, dt = 100

NL Step =  0, fnorm =  7.90072E-03

NL Step =  1, fnorm =  2.01919E-04

NL Step =  2, fnorm =  1.06960E-05

NL Step =  3, fnorm =  2.41683E-09

 Solve Converged!

Notice the obvious ‘stagnant’ in residual for the 2nd order method while not in 
the 1st order.
For the same problem, the wall time is 10 sec vs 6 sec. I would be happy if I 
can reduce 2 sec for the 2nd order method.

-Ling

From: Barry Smith 
Date: Sunday, March 3, 2024 at 12:06 PM
To: Zou, Ling 
Cc: petsc-users@mcs.anl.gov 
Subject: Re: [petsc-users] 'Preconditioning' with lower-order method
Are you forming the Jacobian for the first and second order cases inside of 
Newton? You can run both with -log_view to see how much time is spent in the 
various events (compute function, compute Jacobian, linear solve, .. . ) for 
the two cases
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

   Are you forming the Jacobian for the first and second order cases inside of 
Newton?

   You can run both with -log_view to see how much time is spent in the various 
events (compute function, compute Jacobian, linear solve, ...) for the two 
cases and compare them.




On Mar 3, 2024, at 11:42 AM, Zou, Ling via petsc-users 
 wrote:

Original email may have been sent to the incorrect place.
See below.

-Ling

From: Zou, Ling mailto:l...@anl.gov>>
Date: Sunday, March 3, 2024 at 10:34 AM
To: petsc-users 
mailto:petsc-users-boun...@mcs.anl.gov>>
Subject: 'Preconditioning' with lower-order method
Hi all,

I am solving a PDE system over a spatial domain. Numerical methods are:

  *   Finite volume method (both 1st and 2nd order implemented)
  *   BDF1 and BDF2 for time integration.
What I have noticed is that 1st order FVM converges much faster than 2nd order 
FVM, regardless the time integration scheme. Well, not surprising since 2nd 
order FVM introduces additional non-linearity.

I’m thinking about two possible ways to speed up 2nd order FVM, and would like 
to get some thoughts or community knowledge before jumping into code 
implementation.

Say, let the 2nd order FVM residual function be F2(x) = 0; and the 1st order 
FVM residual function be F1(x) = 0.

  1.  Option – 1, multi-step for each time step
Step 1: solving F1(x) = 0 to obtain a temporary solution x1
Step 2: feed x1 as an initial guess to solve F2(x) = 0 to obtain the final 
solution.
[Not sure if gain any saving at all]


  1.  Option -2, dynamically changing residual function F(x)
In pseudo code, would be something like.

snesFormFunction(SNES snes, Vec u, Vec f, void *)
{
  if (snes.nl_it_no < 4) // 4 being arbitrary here
f = F1(u);
  else
f = F2(u);
}

I know this might be a bit crazy since it may crash after switching residual 
function, still, any thoughts?

Best,

-Ling

Re: [petsc-users] FW: 'Preconditioning' with lower-order method

[Zou, Ling via petsc-users]

Thank you, Mark. This is encouraging!
I will give it a try and report back.

-Ling

From: Mark Adams 
Date: Sunday, March 3, 2024 at 11:10 AM
To: Zou, Ling 
Cc: petsc-users@mcs.anl.gov 
Subject: Re: [petsc-users] FW: 'Preconditioning' with lower-order method
On Sun, Mar 3, 2024 at 11: 42 AM Zou, Ling via petsc-users  wrote: Original email may have been sent to the incorrect place. See 
below. -Ling From: Zou, Ling  Date: Sunday, March 3, 2024 at
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd


On Sun, Mar 3, 2024 at 11:42 AM Zou, Ling via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:
Original email may have been sent to the incorrect place.
See below.

-Ling

From: Zou, Ling mailto:l...@anl.gov>>
Date: Sunday, March 3, 2024 at 10:34 AM
To: petsc-users 
mailto:petsc-users-boun...@mcs.anl.gov>>
Subject: 'Preconditioning' with lower-order method
Hi all,

I am solving a PDE system over a spatial domain. Numerical methods are:

  *   Finite volume method (both 1st and 2nd order implemented)
  *   BDF1 and BDF2 for time integration.
What I have noticed is that 1st order FVM converges much faster than 2nd order 
FVM, regardless the time integration scheme. Well, not surprising since 2nd 
order FVM introduces additional non-linearity.

I’m thinking about two possible ways to speed up 2nd order FVM, and would like 
to get some thoughts or community knowledge before jumping into code 
implementation.

Say, let the 2nd order FVM residual function be F2(x) = 0; and the 1st order 
FVM residual function be F1(x) = 0.

  1.  Option – 1, multi-step for each time step
Step 1: solving F1(x) = 0 to obtain a temporary solution x1
Step 2: feed x1 as an initial guess to solve F2(x) = 0 to obtain the final 
solution.
[Not sure if gain any saving at all]


  1.  Option -2, dynamically changing residual function F(x)

In pseudo code, would be something like.


You can try it. I would doubt the linear version (1) would help. This is 
similar to "defect correction" but not the same.

The nonlinear version (2) is something you could try.
I've seen people switch nonlinear solvers like this but not operators.
You could try it.

Mark

snesFormFunction(SNES snes, Vec u, Vec f, void *)

{

  if (snes.nl_it_no < 4) // 4 being arbitrary here

f = F1(u);

  else

f = F2(u);

}



I know this might be a bit crazy since it may crash after switching residual 
function, still, any thoughts?

Best,

-Ling

[petsc-users] FW: 'Preconditioning' with lower-order method

[Zou, Ling via petsc-users]

Original email may have been sent to the incorrect place.
See below.

-Ling

From: Zou, Ling 
Date: Sunday, March 3, 2024 at 10:34 AM
To: petsc-users 
Subject: 'Preconditioning' with lower-order method
Hi all,

I am solving a PDE system over a spatial domain. Numerical methods are:

  *   Finite volume method (both 1st and 2nd order implemented)
  *   BDF1 and BDF2 for time integration.
What I have noticed is that 1st order FVM converges much faster than 2nd order 
FVM, regardless the time integration scheme. Well, not surprising since 2nd 
order FVM introduces additional non-linearity.

I’m thinking about two possible ways to speed up 2nd order FVM, and would like 
to get some thoughts or community knowledge before jumping into code 
implementation.

Say, let the 2nd order FVM residual function be F2(x) = 0; and the 1st order 
FVM residual function be F1(x) = 0.

  1.  Option – 1, multi-step for each time step
Step 1: solving F1(x) = 0 to obtain a temporary solution x1
Step 2: feed x1 as an initial guess to solve F2(x) = 0 to obtain the final 
solution.
[Not sure if gain any saving at all]


  1.  Option -2, dynamically changing residual function F(x)

In pseudo code, would be something like.



snesFormFunction(SNES snes, Vec u, Vec f, void *)

{

  if (snes.nl_it_no < 4) // 4 being arbitrary here

f = F1(u);

  else

f = F2(u);

}



I know this might be a bit crazy since it may crash after switching residual 
function, still, any thoughts?

Best,

-Ling

Re: [petsc-users] A problem with MatFDColoringSetFunction

[Zou, Ling via petsc-users]

A follow up:

The compiler option '-Wno-bad-function-cast' did not work for c++. See message 
below:
cc1plus: error: command line option '-Wno-bad-function-cast' is valid for 
C/ObjC but not for C++ [-Werror]

Credit to David Andrs of INL, here is an alternative solution:

MatFDColoringSetFunction(fdcoloring,
  (PetscErrorCode(*)(void))(void(*)(void))SNESFormFunction, this);

Thanks,

-Ling

From: petsc-users  on behalf of Zou, Ling via 
petsc-users 
Date: Monday, January 10, 2022 at 9:23 PM
To: Barry Smith 
Cc: PETSc 
Subject: Re: [petsc-users] A problem with MatFDColoringSetFunction
Thank you Barry.

-Ling

From: Barry Smith 
Date: Monday, January 10, 2022 at 5:13 PM
To: Zou, Ling 
Cc: PETSc 
Subject: Re: [petsc-users] A problem with MatFDColoringSetFunction

  This is annoying. You need to pass in a compiler flag to turn off the error 
conditioner for casting a function pointer. It may be 
-Wnoerror=cast-function-type just google.

  Barry


On Jan 10, 2022, at 4:35 PM, Zou, Ling via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:

Hi All,

I would appreciate if you could give some advice for setting the ‘FormFunction’ 
for the MatFDColoringSetFunction function call.
I follow what is shown in the example:
https://petsc.org/release/src/snes/tutorials/ex14.c.html

I setup similar code structure like:
PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
Then use it as:
MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))SNESFormFunction, 
this);

This works fine on my local MacOS. However, when commit it to a remote repo, 
the compiler there gives me the following error (where warnings are treated as 
errors):

src/base/PETScProblemInterface.C:100:67: error:
cast between incompatible function types from
'PetscErrorCode (*)(SNES, Vec, Vec, void*)' {aka 'int (*)(_p_SNES*, _p_Vec*, 
_p_Vec*, void*)'}
to
'PetscErrorCode (*)()' {aka 'int (*)()'} [-Werror=cast-function-type]

Any help is appreciated.

-Ling

Re: [petsc-users] A problem with MatFDColoringSetFunction

[Zou, Ling via petsc-users]

Thank you Barry.

-Ling

From: Barry Smith 
Date: Monday, January 10, 2022 at 5:13 PM
To: Zou, Ling 
Cc: PETSc 
Subject: Re: [petsc-users] A problem with MatFDColoringSetFunction

  This is annoying. You need to pass in a compiler flag to turn off the error 
conditioner for casting a function pointer. It may be 
-Wnoerror=cast-function-type just google.

  Barry



On Jan 10, 2022, at 4:35 PM, Zou, Ling via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:

Hi All,

I would appreciate if you could give some advice for setting the ‘FormFunction’ 
for the MatFDColoringSetFunction function call.
I follow what is shown in the example:
https://petsc.org/release/src/snes/tutorials/ex14.c.html

I setup similar code structure like:
PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
Then use it as:
MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))SNESFormFunction, 
this);

This works fine on my local MacOS. However, when commit it to a remote repo, 
the compiler there gives me the following error (where warnings are treated as 
errors):

src/base/PETScProblemInterface.C:100:67: error:
cast between incompatible function types from
'PetscErrorCode (*)(SNES, Vec, Vec, void*)' {aka 'int (*)(_p_SNES*, _p_Vec*, 
_p_Vec*, void*)'}
to
'PetscErrorCode (*)()' {aka 'int (*)()'} [-Werror=cast-function-type]

Any help is appreciated.

-Ling

[petsc-users] A problem with MatFDColoringSetFunction

[Zou, Ling via petsc-users]

Hi All,

I would appreciate if you could give some advice for setting the ‘FormFunction’ 
for the MatFDColoringSetFunction function call.
I follow what is shown in the example:
https://petsc.org/release/src/snes/tutorials/ex14.c.html

I setup similar code structure like:
PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
Then use it as:
MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))SNESFormFunction, 
this);

This works fine on my local MacOS. However, when commit it to a remote repo, 
the compiler there gives me the following error (where warnings are treated as 
errors):

src/base/PETScProblemInterface.C:100:67: error:
cast between incompatible function types from
'PetscErrorCode (*)(SNES, Vec, Vec, void*)' {aka 'int (*)(_p_SNES*, _p_Vec*, 
_p_Vec*, void*)'}
to
'PetscErrorCode (*)()' {aka 'int (*)()'} [-Werror=cast-function-type]

Any help is appreciated.

-Ling

Re: [petsc-users] Install PETSc 3.14 with Conda type environment

[Zou, Ling via petsc-users]

Barry, you were right. That was the gcc version issue.
This has been resolved with later version of gcc.

-Ling

From: petsc-users 
Date: Sunday, November 15, 2020 at 5:03 PM
To: Barry Smith 
Cc: petsc-users@mcs.anl.gov 
Subject: Re: [petsc-users] Install PETSc 3.14 with Conda type environment
LOL, you are definitely right about the gcc version.
I am still new to conda, and pulled out whatever comes first (through some 
quick google search).

I will update my gcc, and let you know.

-Ling

From: Barry Smith 
Date: Sunday, November 15, 2020 at 4:46 PM
To: Zou, Ling 
Cc: petsc-users@mcs.anl.gov 
Subject: Re: [petsc-users] Install PETSc 3.14 with Conda type environment
 Ling,

 Configure  is failing here.

TESTING: checkCPreprocessor from 
config.setCompilers(config/BuildSystem/config/setCompilers.py:699)
  Locate a functional C preprocessor
Checking for program /Users/lzou/miniconda3/bin/gcc...found
  Defined make macro "CPP" to "gcc -E"
Preprocessing source:
#include "confdefs.h"
#include "conffix.h"
#include 
Executing: gcc -E  
-I/var/folders/fd/yx1jblsj2rzb981kx8ksm03ssqvm65/T/petsc-2jvhontx/config.setCompilers
  
/var/folders/fd/yx1jblsj2rzb981kx8ksm03ssqvm65/T/petsc-2jvhontx/config.setCompilers/conftest.c
Possible ERROR while running preprocessor: exit code 1
stdout:
# 1 
"/var/folders/fd/yx1jblsj2rzb981kx8ksm03ssqvm65/T/petsc-2jvhontx/config.setCompilers/conftest.c"

So for some reason the gcc -E seems to work correctly but produces an non-zero 
exit code?

I notice something very concerning in your Conda environment. It seems to have 
a gcc (GCC) 4.8.5 which is more ancient than my son who is in college. I cannot 
imagine you would want to use such an old compiler. Perhaps if you upgrade the 
Conda environment to use a recent gcc then the other problem will go away.

Please give that a try and let us know how it goes, send configure.log again if 
configure still fails.

  Good luck,

  Barry


  Locate a functional C compiler
Checking for program /Users/lzou/miniconda3/bin/gcc...found
Executing: gcc --version
stdout:
gcc (GCC) 4.8.5
Copyright (C) 2015 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.




On Nov 15, 2020, at 4:27 PM, Zou, Ling via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:

Hi,

I have installed PETSc many times before, so not so new to it.
Recently, I have switched to conda system to manage my working environment, and 
would like to reinstall PETSc, and had some issue during the ‘config’ stage:

***
 UNABLE to CONFIGURE with GIVEN OPTIONS(see configure.log for 
details):
---
Cannot find a C preprocessor
***

Config file is also attached.
Appreciate your help.

-Ling

Re: [petsc-users] Install PETSc 3.14 with Conda type environment

[Zou, Ling via petsc-users]

LOL, you are definitely right about the gcc version.
I am still new to conda, and pulled out whatever comes first (through some 
quick google search).

I will update my gcc, and let you know.

-Ling

From: Barry Smith 
Date: Sunday, November 15, 2020 at 4:46 PM
To: Zou, Ling 
Cc: petsc-users@mcs.anl.gov 
Subject: Re: [petsc-users] Install PETSc 3.14 with Conda type environment
 Ling,

 Configure  is failing here.

TESTING: checkCPreprocessor from 
config.setCompilers(config/BuildSystem/config/setCompilers.py:699)
  Locate a functional C preprocessor
Checking for program /Users/lzou/miniconda3/bin/gcc...found
  Defined make macro "CPP" to "gcc -E"
Preprocessing source:
#include "confdefs.h"
#include "conffix.h"
#include 
Executing: gcc -E  
-I/var/folders/fd/yx1jblsj2rzb981kx8ksm03ssqvm65/T/petsc-2jvhontx/config.setCompilers
  
/var/folders/fd/yx1jblsj2rzb981kx8ksm03ssqvm65/T/petsc-2jvhontx/config.setCompilers/conftest.c
Possible ERROR while running preprocessor: exit code 1
stdout:
# 1 
"/var/folders/fd/yx1jblsj2rzb981kx8ksm03ssqvm65/T/petsc-2jvhontx/config.setCompilers/conftest.c"

So for some reason the gcc -E seems to work correctly but produces an non-zero 
exit code?

I notice something very concerning in your Conda environment. It seems to have 
a gcc (GCC) 4.8.5 which is more ancient than my son who is in college. I cannot 
imagine you would want to use such an old compiler. Perhaps if you upgrade the 
Conda environment to use a recent gcc then the other problem will go away.

Please give that a try and let us know how it goes, send configure.log again if 
configure still fails.

  Good luck,

  Barry


  Locate a functional C compiler
Checking for program /Users/lzou/miniconda3/bin/gcc...found
Executing: gcc --version
stdout:
gcc (GCC) 4.8.5
Copyright (C) 2015 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.



On Nov 15, 2020, at 4:27 PM, Zou, Ling via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:

Hi,

I have installed PETSc many times before, so not so new to it.
Recently, I have switched to conda system to manage my working environment, and 
would like to reinstall PETSc, and had some issue during the ‘config’ stage:

***
 UNABLE to CONFIGURE with GIVEN OPTIONS(see configure.log for 
details):
---
Cannot find a C preprocessor
***

Config file is also attached.
Appreciate your help.

-Ling

Re: [petsc-users] [EXT]Re: Discontinuities in the Jacobian matrix for nonlinear problem

[Zou, Ling via petsc-users]

Ok, I do not pretend to be shallow water expert here. I only played with it 
maybe ten years ago.

When you say 5 point stencil in 1d flow, that sounds so much like a finite 
volume method with TVD reconstruction.
I solved many 1d flow equation using PETSc (Euler eqn., single-phase flow eqn., 
two-phase flow eqn.), had a lot of ‘zero-flow’ conditions, so far PETSc works 
pretty well for me.

-Ling

From: Alexander B Prescott 
Date: Monday, April 6, 2020 at 8:44 PM
To: Nathan Collier 
Cc: "Zou, Ling" , PETSc 
Subject: Re: [EXT]Re: [petsc-users] Discontinuities in the Jacobian matrix for 
nonlinear problem

Hi Nathan,

Thanks for the thoughtful response. For the 1D toy version I force flow to 
occur at every node, but my ultimate objective is to apply this to a 2D 
floodplain in which a node may or may not have flow in the 'true' solution. 
Based on what you point out, however, it seems that this approach may be doomed 
to fail. For one thing, using a 5 point stencil, a node with no flow and whose 
neighbor's have no flow will produce a row in the 'bent physics' Jacobian of 
all zero entries. Based on my experience so far that will cause PETSc to return 
errors. Does that sound right to you?

I agree that the difficulty encountered with this sort of problem is 
counterintuitive to the relative simplification made on nature!

Best,
Alexander

On Mon, Apr 6, 2020 at 12:45 PM Nathan Collier 
mailto:nathaniel.coll...@gmail.com>> wrote:

External Email
Alexander,

I am not familiar with your specific model, but I do have experience working 
with the diffusive / kinematic wave approximation to shallow water equations. I 
always had a lot of trouble near ponding conditions (steady state), which is 
odd because when nature says do nothing, you don't expect the equations to be 
hard to solve. Some kind soul pointed out to me that this was because the 
equations are derived using a power law model relating flow to slope/depth 
which breaks down when you should have no flow. So as long as you have flow, 
the solver behaves well and things are fine (i.e. smooth terrain, academic test 
problems), but when you have local ponding due to real terrain effects, you are 
hosed. A better solution approach doesn't help if your equations aren't meant 
for the problem you are trying to solve.

I think 'bending' the physics is the right idea, but I never had much luck 
myself with this.

Nate






On Mon, Apr 6, 2020 at 2:56 PM Zou, Ling via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:
What about ‘bending’ the physics a bit?

Q_i-1/2 proportional to sqrt(x[i-1] + z[i-1] - (x[i] + z[i])),  
if x[i-1]+z[i-1]  >  x[i]+z[i] + eps
Q_i-1/2 proportional to -1.0*sqrt(x[i] + z[i] - (x[i-1] + z[i-1])), if  
  x[i]+z[i]  >  x[i-1]+z[i-1] + eps
Q_i-1/2 proportional to x[i-1] + z[i-1] - (x[i] + z[i]) 
  in between

in which, eps is a very small positive number.

-Ling


From: petsc-users 
mailto:petsc-users-boun...@mcs.anl.gov>> on 
behalf of Alexander B Prescott 
mailto:alexpresc...@email.arizona.edu>>
Date: Monday, April 6, 2020 at 1:06 PM
To: PETSc mailto:petsc-users@mcs.anl.gov>>
Subject: [petsc-users] Discontinuities in the Jacobian matrix for nonlinear 
problem

Hello,

The non-linear boundary-value problem I am applying PETSc to is a relatively 
simple steady-state flow routing algorithm based on the continuity equation, 
such that Div(Q) = 0 everywhere (Q=discharge). I use a finite volume approach 
to calculate flow between nodes, with Q calculated as a piecewise smooth 
function of the local flow depth and the water-surface slope. In 1D, the 
residual is calculated as R(x_i)=Q_i-1/2 - Q_i+1/2.
For example, Q_i-1/2 at x[i]:
Q_i-1/2 proportional to sqrt(x[i-1] + z[i-1] - (x[i] + z[i])),  if  
x[i-1]+z[i-1]  >  x[i]+z[i]
Q_i-1/2 proportional to -1.0*sqrt(x[i] + z[i] - (x[i-1] + z[i-1])),  if 
x[i]+z[i]  >  x[i-1]+z[i-1]

Where z[i] is local topography and doesn't change over the iterations, and 
Q_i+1/2 is computed analogously. So the residual derivatives with respect to 
x[i-1], x[i] and x[i+1] are not continuous when the water-surface slope = 0.

Are there intelligent ways to handle this problem? My 1D trial runs naively fix 
any zero-valued water-surface slopes to a small non-zero positive value (e.g. 
1e-12). Solver convergence has been mixed and highly dependent on the initial 
guess. So far, FAS with QN coarse solver has been the most robust.

Restricting x[i] to be non-negative is a separate issue, to which I have 
applied the SNES_VI solvers. They perform modestly but have been less robust.

Best,
Alexander



--
Alexander Prescott
alexpresc...@email.arizona.edu<mailto:alexpresc...@email.arizona.edu>
PhD Candidate, The University of Arizona
Department of Geosciences
1040 E. 4th Street
Tucson, AZ, 85721


--
Alexander Prescott
alexpresc...@email.arizona.edu<mailto:

Re: [petsc-users] Discontinuities in the Jacobian matrix for nonlinear problem

[Zou, Ling via petsc-users]

What about ‘bending’ the physics a bit?

Q_i-1/2 proportional to sqrt(x[i-1] + z[i-1] - (x[i] + z[i])),  
if x[i-1]+z[i-1]  >  x[i]+z[i] + eps
Q_i-1/2 proportional to -1.0*sqrt(x[i] + z[i] - (x[i-1] + z[i-1])), if  
  x[i]+z[i]  >  x[i-1]+z[i-1] + eps
Q_i-1/2 proportional to x[i-1] + z[i-1] - (x[i] + z[i]) 
  in between

in which, eps is a very small positive number.

-Ling


From: petsc-users  on behalf of Alexander B 
Prescott 
Date: Monday, April 6, 2020 at 1:06 PM
To: PETSc 
Subject: [petsc-users] Discontinuities in the Jacobian matrix for nonlinear 
problem

Hello,

The non-linear boundary-value problem I am applying PETSc to is a relatively 
simple steady-state flow routing algorithm based on the continuity equation, 
such that Div(Q) = 0 everywhere (Q=discharge). I use a finite volume approach 
to calculate flow between nodes, with Q calculated as a piecewise smooth 
function of the local flow depth and the water-surface slope. In 1D, the 
residual is calculated as R(x_i)=Q_i-1/2 - Q_i+1/2.
For example, Q_i-1/2 at x[i]:
Q_i-1/2 proportional to sqrt(x[i-1] + z[i-1] - (x[i] + z[i])),  if  
x[i-1]+z[i-1]  >  x[i]+z[i]
Q_i-1/2 proportional to -1.0*sqrt(x[i] + z[i] - (x[i-1] + z[i-1])),  if 
x[i]+z[i]  >  x[i-1]+z[i-1]

Where z[i] is local topography and doesn't change over the iterations, and 
Q_i+1/2 is computed analogously. So the residual derivatives with respect to 
x[i-1], x[i] and x[i+1] are not continuous when the water-surface slope = 0.

Are there intelligent ways to handle this problem? My 1D trial runs naively fix 
any zero-valued water-surface slopes to a small non-zero positive value (e.g. 
1e-12). Solver convergence has been mixed and highly dependent on the initial 
guess. So far, FAS with QN coarse solver has been the most robust.

Restricting x[i] to be non-negative is a separate issue, to which I have 
applied the SNES_VI solvers. They perform modestly but have been less robust.

Best,
Alexander



--
Alexander Prescott
alexpresc...@email.arizona.edu
PhD Candidate, The University of Arizona
Department of Geosciences
1040 E. 4th Street
Tucson, AZ, 85721

Re: [petsc-users] Problem about Residual evaluation

[Zou, Ling via petsc-users]

For question 2, there are some discussions in PETSc’s manual section 4.3 and 
4.3.3.
https://www.mcs.anl.gov/petsc/petsc-current/docs/manual.pdf

-Ling

From: petsc-users  on behalf of Yingjie Wu via 
petsc-users 
Reply-To: Yingjie Wu 
Date: Monday, April 1, 2019 at 9:22 AM
To: "petsc-users@mcs.anl.gov" 
Subject: [petsc-users] Problem about Residual evaluation

Dear PETSc developers:
Hi,

I've been using -snes_mf_operator and I've customized a precondition matrix to 
solve my problem.I have two questions about the residuals of linear steps(KSP 
residual).



1.Since I'm using a matrix-free method, how do we get KSP residuals in PETSc?

r_m = b - A*x_m

Is finite difference used to approximate "A*x_m" ?



2.What is the difference between instruction ' -ksp_monitor ' and ' 
-ksp_monitor_true_residual ' in how they are calculated?



Thanks,

Yingjie