[petsc-users] Diagnosing Poisson Solver Behavior

[K. N. Ramachandran]

Hello all,

I am a graduate student pursuing my Master's and I am trying to benchmark a
previous work by using PETSc for solving Poisson's Equation.

I am starting off with a serial code and I am trying to keep my code
modular, i.e. I generate the sparse matrix format and send it to PETSc or
any other solver. So I haven't built my code from the ground up using
PETSc's native data structures.

I am having trouble understanding the behavior of the solver and would like
your thoughts or inputs on what I can do better. I have both Dirichlet and
Neumann boundary conditions and my matrix size (number of rows) is around a
million but very sparse (~7 nonzeros per row), as can be expected from a
finite difference discretization of Poisson's equation.

I tried the methods outlined here
<http://scicomp.stackexchange.com/questions/513/why-is-my-iterative-linear-solver-not-converging?rq=1>
and here
<http://scicomp.stackexchange.com/questions/34/how-can-i-estimate-the-condition-number-of-a-large-sparse-matrix-using-petsc>.
Reverting to a 41^3 grid, I got the approximate condition number
(using -ksp_monitor_singular_value
-ksp_type gmres -ksp_gmres_restart 1000 -pc_type none) as ~9072, which
seems pretty large. Higher matrix sizes give a larger condition number.

1) My best performing solver + preconditioner is bcgs+ilu(0) (on 1e6 grid)
which solves in around 32 seconds, 196 iterations. How do I get a fix for
what the lower bound on the running time could be?

2) Initially -pc_type hypre just Diverged and I was never able to use it.
Looking at this thread
<http://lists.mcs.anl.gov/pipermail/petsc-users/2013-October/019127.html>,
I had tried the options and it no longer diverges, but the residuals reduce
and hover around a constant value. How do I work with hypre to get a useful
preconditioner?

Initially I solve Laplace's equation, so the mesh grid size has no effect
and even when I solve Poisson's equation, the spacing is carried over to
the RHS, so I am pretty sure the spacing is not affecting the condition
number calculation.

Hope this helps. Please let me know if you might need more information.

Thanking You,
K.N.Ramachandran
Ph: 814-441-4279

[K. N. Ramachandran]

Sorry, some more questions.

3) Also, for Dirichlet bc, I specify the value through Identity rows, i.e.
A_ii = 1 and the rhs value would correspond to the Dirichlet condition. I
am specifying it this way for my convenience. I am aware that
MatZerosRowColumns might help here, but would keeping it this way be
detrimental?

4) Can I expect symmetric matrices to perform better, i.e. if I eliminate
Dirichlet rows? But I would still be left with Neumann boundary conditions,
where I use the second order formulation. If I used the first order
formulation and made it symmetric, would that be an advantage? I tried the
latter, but I didn't see the condition number change much.

Post by K. N. Ramachandran
Hello all,
I am a graduate student pursuing my Master's and I am trying to benchmark
a previous work by using PETSc for solving Poisson's Equation.
I am starting off with a serial code and I am trying to keep my code
modular, i.e. I generate the sparse matrix format and send it to PETSc or
any other solver. So I haven't built my code from the ground up using
PETSc's native data structures.
I am having trouble understanding the behavior of the solver and would
like your thoughts or inputs on what I can do better. I have both Dirichlet
and Neumann boundary conditions and my matrix size (number of rows) is
around a million but very sparse (~7 nonzeros per row), as can be expected
from a finite difference discretization of Poisson's equation.
I tried the methods outlined here
<http://scicomp.stackexchange.com/questions/513/why-is-my-iterative-linear-solver-not-converging?rq=1>
and here
<http://scicomp.stackexchange.com/questions/34/how-can-i-estimate-the-condition-number-of-a-large-sparse-matrix-using-petsc>.
Reverting to a 41^3 grid, I got the approximate condition number (using -ksp_monitor_singular_value
-ksp_type gmres -ksp_gmres_restart 1000 -pc_type none) as ~9072, which
seems pretty large. Higher matrix sizes give a larger condition number.
1) My best performing solver + preconditioner is bcgs+ilu(0) (on 1e6
grid) which solves in around 32 seconds, 196 iterations. How do I get a
fix for what the lower bound on the running time could be?
2) Initially -pc_type hypre just Diverged and I was never able to use it.
Looking at this thread
<http://lists.mcs.anl.gov/pipermail/petsc-users/2013-October/019127.html>,
I had tried the options and it no longer diverges, but the residuals reduce
and hover around a constant value. How do I work with hypre to get a
useful preconditioner?
Initially I solve Laplace's equation, so the mesh grid size has no effect
and even when I solve Poisson's equation, the spacing is carried over to
the RHS, so I am pretty sure the spacing is not affecting the condition
number calculation.
Hope this helps. Please let me know if you might need more information.
Thanking You,
K.N.Ramachandran
Ph: 814-441-4279

Thanking You,
K.N.Ramachandran
Ph: 814-441-4279

[Matthew Knepley]

Post by K. N. Ramachandran
Sorry, some more questions.
3) Also, for Dirichlet bc, I specify the value through Identity rows, i.e.
A_ii = 1 and the rhs value would correspond to the Dirichlet condition. I
am specifying it this way for my convenience. I am aware that
MatZerosRowColumns might help here, but would keeping it this way be
detrimental?

That is fine. However you would want to scale these entries to be
approximately the same size as the other diagonal entries.

Post by K. N. Ramachandran
4) Can I expect symmetric matrices to perform better, i.e. if I eliminate
Dirichlet rows? But I would still be left with Neumann boundary conditions,
where I use the second order formulation. If I used the first order
formulation and made it symmetric, would that be an advantage? I tried the
latter, but I didn't see the condition number change much.

Will not matter for MG.

Matt

Post by K. N. Ramachandran

Post by K. N. Ramachandran
Hello all,
I am a graduate student pursuing my Master's and I am trying to benchmark
a previous work by using PETSc for solving Poisson's Equation.
I am starting off with a serial code and I am trying to keep my code
modular, i.e. I generate the sparse matrix format and send it to PETSc or
any other solver. So I haven't built my code from the ground up using
PETSc's native data structures.
I am having trouble understanding the behavior of the solver and would
like your thoughts or inputs on what I can do better. I have both Dirichlet
and Neumann boundary conditions and my matrix size (number of rows) is
around a million but very sparse (~7 nonzeros per row), as can be expected
from a finite difference discretization of Poisson's equation.
I tried the methods outlined here
<http://scicomp.stackexchange.com/questions/513/why-is-my-iterative-linear-solver-not-converging?rq=1>
and here
<http://scicomp.stackexchange.com/questions/34/how-can-i-estimate-the-condition-number-of-a-large-sparse-matrix-using-petsc>.
Reverting to a 41^3 grid, I got the approximate condition number (using -ksp_monitor_singular_value
-ksp_type gmres -ksp_gmres_restart 1000 -pc_type none) as ~9072, which
seems pretty large. Higher matrix sizes give a larger condition number.
1) My best performing solver + preconditioner is bcgs+ilu(0) (on 1e6
grid) which solves in around 32 seconds, 196 iterations. How do I get a
fix for what the lower bound on the running time could be?
2) Initially -pc_type hypre just Diverged and I was never able to use
it. Looking at this thread
<http://lists.mcs.anl.gov/pipermail/petsc-users/2013-October/019127.html>,
I had tried the options and it no longer diverges, but the residuals reduce
and hover around a constant value. How do I work with hypre to get a
useful preconditioner?
Initially I solve Laplace's equation, so the mesh grid size has no effect
and even when I solve Poisson's equation, the spacing is carried over to
the RHS, so I am pretty sure the spacing is not affecting the condition
number calculation.
Hope this helps. Please let me know if you might need more information.
Thanking You,
K.N.Ramachandran
Ph: 814-441-4279

Thanking You,
K.N.Ramachandran
Ph: 814-441-4279

[K. N. Ramachandran]

Hello Matt,

Thanks for your reply. I am actually retrying against the latest PETSc and
will revert once I get the information.

Post by Matthew Knepley

Post by K. N. Ramachandran
Sorry, some more questions.
3) Also, for Dirichlet bc, I specify the value through Identity rows,
i.e. A_ii = 1 and the rhs value would correspond to the Dirichlet
condition. I am specifying it this way for my convenience. I am aware that
MatZerosRowColumns might help here, but would keeping it this way be
detrimental?

That is fine. However you would want to scale these entries to be
approximately the same size as the other diagonal entries.

Post by K. N. Ramachandran
4) Can I expect symmetric matrices to perform better, i.e. if I eliminate
Dirichlet rows? But I would still be left with Neumann boundary conditions,
where I use the second order formulation. If I used the first order
formulation and made it symmetric, would that be an advantage? I tried the
latter, but I didn't see the condition number change much.

Will not matter for MG.
Matt

Regards,
K.N.Ramachandran
Ph: 814-441-4279

[K. N. Ramachandran]

Hello Matt,

Actually I felt the boundary conditions were having a role to play and set
all the boundary conditions to Dirichlet. In this case, convergence was
almost immediate with the Hypre preconditioner, taking 17 seconds with 3
iterations. The MG method took about the same time though.

So I reverted to the Dirichlet, Neumann mix of BCs and Hypre starts to
diverge. Please find attached the output for the Hypre run using Dirichlet
and Neumann for a 21^3 grid (rows), with a max of 7 nonzeros per row.
Details of the options used before running are in the file. The solver used
in all cases is bcgs.

Also attached is the MG output for 101^3 grid
1) Dirichlet and Neumann
2) Dirichlet only

where it seems to take about the same time.

I was thinking if the Null space has something to do with this? From the
PETSc docs, I understand that the Null space should be removed/attached
when solving singular systems. My domain will have at least two Dirichlet
conditions, so it is not singular. But since the Neumann conditions seem to
affect the convergence with Hypre, perhaps they have a role to play in
terms of the null space?

Post by K. N. Ramachandran
Hello Matt,
Thanks for your reply. I am actually retrying against the latest PETSc and
will revert once I get the information.

Post by Matthew Knepley

Post by K. N. Ramachandran
Sorry, some more questions.
3) Also, for Dirichlet bc, I specify the value through Identity rows,
i.e. A_ii = 1 and the rhs value would correspond to the Dirichlet
condition. I am specifying it this way for my convenience. I am aware that
MatZerosRowColumns might help here, but would keeping it this way be
detrimental?

That is fine. However you would want to scale these entries to be
approximately the same size as the other diagonal entries.

Post by K. N. Ramachandran
4) Can I expect symmetric matrices to perform better, i.e. if I
eliminate Dirichlet rows? But I would still be left with Neumann boundary
conditions, where I use the second order formulation. If I used the first
order formulation and made it symmetric, would that be an advantage? I
tried the latter, but I didn't see the condition number change much.

Will not matter for MG.
Matt

Regards,
K.N.Ramachandran
Ph: 814-441-4279

Thanking You,
K.N.Ramachandran
Ph: 814-441-4279

[Matthew Knepley]

Post by K. N. Ramachandran
Hello Matt,
Actually I felt the boundary conditions were having a role to play and set
all the boundary conditions to Dirichlet. In this case, convergence was
almost immediate with the Hypre preconditioner, taking 17 seconds with 3
iterations. The MG method took about the same time though.
So I reverted to the Dirichlet, Neumann mix of BCs and Hypre starts to
diverge. Please find attached the output for the Hypre run using Dirichlet
and Neumann for a 21^3 grid (rows), with a max of 7 nonzeros per row.
Details of the options used before running are in the file. The solver used
in all cases is bcgs.
Also attached is the MG output for 101^3 grid
1) Dirichlet and Neumann
2) Dirichlet only
where it seems to take about the same time.

Notice that you have no levels of MG here. You need to use -pc_mg_levels <n>

Matt

Post by K. N. Ramachandran
I was thinking if the Null space has something to do with this? From the
PETSc docs, I understand that the Null space should be removed/attached
when solving singular systems. My domain will have at least two Dirichlet
conditions, so it is not singular. But since the Neumann conditions seem to
affect the convergence with Hypre, perhaps they have a role to play in
terms of the null space?

Post by K. N. Ramachandran
Hello Matt,
Thanks for your reply. I am actually retrying against the latest PETSc
and will revert once I get the information.

Post by Matthew Knepley

Post by K. N. Ramachandran
Sorry, some more questions.
3) Also, for Dirichlet bc, I specify the value through Identity rows,
i.e. A_ii = 1 and the rhs value would correspond to the Dirichlet
condition. I am specifying it this way for my convenience. I am aware that
MatZerosRowColumns might help here, but would keeping it this way be
detrimental?

That is fine. However you would want to scale these entries to be
approximately the same size as the other diagonal entries.

Post by K. N. Ramachandran
4) Can I expect symmetric matrices to perform better, i.e. if I
eliminate Dirichlet rows? But I would still be left with Neumann boundary
conditions, where I use the second order formulation. If I used the first
order formulation and made it symmetric, would that be an advantage? I
tried the latter, but I didn't see the condition number change much.

Will not matter for MG.
Matt

Regards,
K.N.Ramachandran
Ph: 814-441-4279

Thanking You,
K.N.Ramachandran
Ph: 814-441-4279

[K. N. Ramachandran]

Hello Matt,

Ah ok. I had tried specifying -pc_mg_levels earlier, but it asks for
interpolation to be setup. I am currently working on using DMDA from
examples 32 and 34, which sets up the interpolation, using
DMDASetInterpolationType. So this is the recommended way to go forward
right?

Since I supply the sparse matrix format and the rhs values to PETSc, does
it make sense to do a DMDACreate1d instead of 3d? I did a sample
implementation using DMDA 1d and when I call KSPSolve and fetch the
solution, it comes out all zeros, with the solver having converged. I am
trying to understand what I am doing wrong there.

Post by Matthew Knepley

Post by K. N. Ramachandran
Hello Matt,
Actually I felt the boundary conditions were having a role to play and
set all the boundary conditions to Dirichlet. In this case, convergence was
almost immediate with the Hypre preconditioner, taking 17 seconds with 3
iterations. The MG method took about the same time though.
So I reverted to the Dirichlet, Neumann mix of BCs and Hypre starts to
diverge. Please find attached the output for the Hypre run using Dirichlet
and Neumann for a 21^3 grid (rows), with a max of 7 nonzeros per row.
Details of the options used before running are in the file. The solver used
in all cases is bcgs.
Also attached is the MG output for 101^3 grid
1) Dirichlet and Neumann
2) Dirichlet only
where it seems to take about the same time.

Notice that you have no levels of MG here. You need to use -pc_mg_levels <n>
Matt

Thanking You,
Ramachandran K.N.
Ph: 814-441-4279

[Matthew Knepley]

Post by K. N. Ramachandran
Hello Matt,
Ah ok. I had tried specifying -pc_mg_levels earlier, but it asks for
interpolation to be setup. I am currently working on using DMDA from
examples 32 and 34, which sets up the interpolation, using
DMDASetInterpolationType. So this is the recommended way to go forward
right?
Since I supply the sparse matrix format and the rhs values to PETSc, does
it make sense to do a DMDACreate1d instead of 3d? I did a sample
implementation using DMDA 1d and when I call KSPSolve and fetch the
solution, it comes out all zeros, with the solver having converged. I am
trying to understand what I am doing wrong there.

My understanding was that you are just solving 3D Poisson with a star
stencil on a regular grid. SNES ex5 solves 2D Poisson with
a star stencil on a regular grid (set lambda = 0.0). Just use that until
you understand the performance, and then change it for 3D
which is only about 10 lines of code.

For this problem, there is no need for the more complicated stuff.

Thanks,

Matt

Post by K. N. Ramachandran

Post by Matthew Knepley

Post by K. N. Ramachandran
Hello Matt,
Actually I felt the boundary conditions were having a role to play and
set all the boundary conditions to Dirichlet. In this case, convergence was
almost immediate with the Hypre preconditioner, taking 17 seconds with 3
iterations. The MG method took about the same time though.
So I reverted to the Dirichlet, Neumann mix of BCs and Hypre starts to
diverge. Please find attached the output for the Hypre run using Dirichlet
and Neumann for a 21^3 grid (rows), with a max of 7 nonzeros per row.
Details of the options used before running are in the file. The solver used
in all cases is bcgs.
Also attached is the MG output for 101^3 grid
1) Dirichlet and Neumann
2) Dirichlet only
where it seems to take about the same time.

Notice that you have no levels of MG here. You need to use -pc_mg_levels <n>
Matt

Thanking You,
Ramachandran K.N.
Ph: 814-441-4279

[Barry Smith]

Post by K. N. Ramachandran
Hello Matt,
Ah ok. I had tried specifying -pc_mg_levels earlier, but it asks for interpolation to be setup. I am currently working on using DMDA from examples 32 and 34, which sets up the interpolation, using DMDASetInterpolationType. So this is the recommended way to go forward right?
Since I supply the sparse matrix format and the rhs values to PETSc, does it make sense to do a DMDACreate1d instead of 3d? I did a sample implementation using DMDA 1d and when I call KSPSolve and fetch the solution, it comes out all zeros, with the solver having converged. I am trying to understand what I am doing wrong there.
My understanding was that you are just solving 3D Poisson with a star stencil on a regular grid. SNES ex5 solves 2D Poisson with
a star stencil on a regular grid (set lambda = 0.0). Just use that until you understand the performance, and then change it for 3D
which is only about 10 lines of code.

See src/ksp/ksp/examples/tutorials/ex45.c which does exactly the 3d problem.

Post by K. N. Ramachandran
For this problem, there is no need for the more complicated stuff.
Thanks,
Matt
Hello Matt,
Actually I felt the boundary conditions were having a role to play and set all the boundary conditions to Dirichlet. In this case, convergence was almost immediate with the Hypre preconditioner, taking 17 seconds with 3 iterations. The MG method took about the same time though.
So I reverted to the Dirichlet, Neumann mix of BCs and Hypre starts to diverge. Please find attached the output for the Hypre run using Dirichlet and Neumann for a 21^3 grid (rows), with a max of 7 nonzeros per row. Details of the options used before running are in the file. The solver used in all cases is bcgs.
Also attached is the MG output for 101^3 grid
1) Dirichlet and Neumann
2) Dirichlet only
where it seems to take about the same time.
Notice that you have no levels of MG here. You need to use -pc_mg_levels <n>
Matt
Thanking You,
Ramachandran K.N.
Ph: 814-441-4279
--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

[K. N. Ramachandran]

Hello all,

@Matt: Actually I can't fully adopt the SNES ex5 since I think I am much
closer to getting the DMDA stuff to work, even if it is just serial. So I
am currently focussing on that.

@Barry: The matrix is coming from the Finite difference discretization of a
3D regular grid. I had a look at KSP ex45.c and I am struggling with what
looks to be a simple problem.

I currently supply the matrix to PETSc in a CSR form. I have specified
ComputeRHS and ComputeMatrix functions, analogous to ex45.c.

If my 3D grid is regular cube with 5 nodes on each side, I create DMDA
using DMDACreate3d((comm world), boundary types, STAR_STENCIL, etc., M = 5,
N = 5, P = 5, dof=1, s=1, PETSC_NULL on remaining).

In the ComputeMatrix function, I try to reuse my CSR form by doing:

#undef __FUNCT__
#define __FUNCT__ "ComputeMatrix"
PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat A, void* ctx)
{
PetscFunctionBeginUser;

Mat temp;
// build up the matrix
MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, numRows, numRows, rows,
cols, values, &temp);

MatCopy(temp, A, DIFFERENT_NONZERO_PATTERN);

MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

PetscFunctionReturn(0);

} // end of ComputeMatrix

Where rows, cols, values are stored globally, for now. Also I am keeping
the code serial, so I am just trying to use DMDA as a convenient way to
take care of interpolation between MG levels, if my understanding is
correct.

But on the MatCopy line, I get an error

DM Object: 1 MPI processes
type: da
Processor [0] M 5 N 5 P 5 m 1 n 1 p 1 w 1 s 1
X range of indices: 0 5, Y range of indices: 0 5, Z range of indices: 0 5
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Inserting a new nonzero at (25,27) in the matrix
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.2, Oct, 02, 2015

I have been trying but I have not been able to resolve this quickly. I feel
like I am creating the DMDA correctly. Is there a better way of reusing a
CSR form to setup a DM Matrix? Should I be doing DMCreateMatrix?

Thanks a lot for your help.


Regards,
K.N.Ramachandran
Ph: 814-441-4279

[Barry Smith]

we need the full error message, what function is causing the error? In theory what you are doing sounds like it should work, it is all the little details that are likely causing the issues. In parallel you will not be able to reuse some old code to "compute the Finite difference discretization of a 3D regular grid" so you might consider just discarding that old code and using what PETSc provides for computing the "Finite difference discretization of a 3D regular grid."

Barry

Post by K. N. Ramachandran
Hello all,
@Matt: Actually I can't fully adopt the SNES ex5 since I think I am much closer to getting the DMDA stuff to work, even if it is just serial. So I am currently focussing on that.
@Barry: The matrix is coming from the Finite difference discretization of a 3D regular grid. I had a look at KSP ex45.c and I am struggling with what looks to be a simple problem.
I currently supply the matrix to PETSc in a CSR form. I have specified ComputeRHS and ComputeMatrix functions, analogous to ex45.c.
If my 3D grid is regular cube with 5 nodes on each side, I create DMDA using DMDACreate3d((comm world), boundary types, STAR_STENCIL, etc., M = 5, N = 5, P = 5, dof=1, s=1, PETSC_NULL on remaining).
#undef __FUNCT__
#define __FUNCT__ "ComputeMatrix"
PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat A, void* ctx)
{
PetscFunctionBeginUser;
Mat temp;
// build up the matrix
MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, numRows, numRows, rows, cols, values, &temp);
MatCopy(temp, A, DIFFERENT_NONZERO_PATTERN);
MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
PetscFunctionReturn(0);
} // end of ComputeMatrix
Where rows, cols, values are stored globally, for now. Also I am keeping the code serial, so I am just trying to use DMDA as a convenient way to take care of interpolation between MG levels, if my understanding is correct.
But on the MatCopy line, I get an error
DM Object: 1 MPI processes
type: da
Processor [0] M 5 N 5 P 5 m 1 n 1 p 1 w 1 s 1
X range of indices: 0 5, Y range of indices: 0 5, Z range of indices: 0 5
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Inserting a new nonzero at (25,27) in the matrix
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.2, Oct, 02, 2015
I have been trying but I have not been able to resolve this quickly. I feel like I am creating the DMDA correctly. Is there a better way of reusing a CSR form to setup a DM Matrix? Should I be doing DMCreateMatrix?
Thanks a lot for your help.
Regards,
K.N.Ramachandran
Ph: 814-441-4279

[K. N. Ramachandran]

Hello Barry,

The full error message is

[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Inserting a new nonzero at (25,27) in the matrix
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.2, Oct, 02, 2015
[0]PETSC ERROR: /gpfs/work/(user name)/progs/PICC_Code/piccode on a
linux-debug named (host) by (user) Thu Oct 15 11:22:41 2015
[0]PETSC ERROR: Configure options --with-fc=mpif90 PETSC_ARCH=linux-debug
--with-debugging=1 --with-clanguage=cxx --with-cxx=mpicxx --with-cc=mpicc
--with-blas-lib="-L/usr/global/intel/mkl/11.0.0.079/lib/intel64 -lmkl_rt"
--with-lapack-lib="-L/usr/global/intel/mkl/11.0.0.079/lib/intel64 -lmkl_rt"
--with-shared-libraries --with-x=0 --download-hypre=1
[0]PETSC ERROR: #1 MatSetValues_SeqAIJ() line 484 in /gpfs/scratch/(user
name)/petsc-3.6.2/src/mat/impls/aij/seq/aij.c
[0]PETSC ERROR: #2 MatSetValues() line 1173 in /gpfs/scratch/(user
name)/petsc-3.6.2/src/mat/interface/matrix.c
[0]PETSC ERROR: #3 MatCopy_Basic() line 3722 in /gpfs/scratch/(user
name)/petsc-3.6.2/src/mat/interface/matrix.c
[0]PETSC ERROR: #4 MatCopy_SeqAIJ() line 2630 in /gpfs/scratch/(user
name)/petsc-3.6.2/src/mat/impls/aij/seq/aij.c
[0]PETSC ERROR: #5 MatCopy() line 3779 in /gpfs/scratch/(user
name)/petsc-3.6.2/src/mat/interface/matrix.c
48 MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
(gdb) n
49 MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);


I tried a DMDACreate1D where I specify the dimension as n^3, where n is the
number of nodes on one side of the cube, but I still get the same error.
The error is generated at the MatCopy line as mentioned earlier.

In parallel you will not be able to reuse some old code to "compute the

Post by Barry Smith
Finite difference discretization of a 3D regular grid" so you might
consider just discarding that old code and using what PETSc provides for
computing the "Finite difference discretization of a 3D regular grid."

Ah ok, I had meant to say Finite difference discretization of the Poisson
equation.

If my parallel code could generate the sparse matrices formats for the
portion owned by each ranks, would there be a convenient way to generate
the PETSc matrices to be part of the DM? I had noticed that
CreateMatSeqAIJWithArrays changes the location of the Mat pointer passed in
for ComputeMatrix, so the actual matrix in DM remains as all zeros, so my
final solution comes out as zeros too.


Thanking You,
K.N.Ramachandran
Ph: 814-441-4279

[Matthew Knepley]

Post by K. N. Ramachandran
Hello Barry,
The full error message is
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Inserting a new nonzero at (25,27) in the matrix
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.2, Oct, 02, 2015
[0]PETSC ERROR: /gpfs/work/(user name)/progs/PICC_Code/piccode on a
linux-debug named (host) by (user) Thu Oct 15 11:22:41 2015
[0]PETSC ERROR: Configure options --with-fc=mpif90 PETSC_ARCH=linux-debug
--with-debugging=1 --with-clanguage=cxx --with-cxx=mpicxx --with-cc=mpicc
--with-blas-lib="-L/usr/global/intel/mkl/11.0.0.079/lib/intel64 -lmkl_rt"
--with-lapack-lib="-L/usr/global/intel/mkl/11.0.0.079/lib/intel64
-lmkl_rt" --with-shared-libraries --with-x=0 --download-hypre=1
[0]PETSC ERROR: #1 MatSetValues_SeqAIJ() line 484 in /gpfs/scratch/(user
name)/petsc-3.6.2/src/mat/impls/aij/seq/aij.c
[0]PETSC ERROR: #2 MatSetValues() line 1173 in /gpfs/scratch/(user
name)/petsc-3.6.2/src/mat/interface/matrix.c
[0]PETSC ERROR: #3 MatCopy_Basic() line 3722 in /gpfs/scratch/(user
name)/petsc-3.6.2/src/mat/interface/matrix.c
[0]PETSC ERROR: #4 MatCopy_SeqAIJ() line 2630 in /gpfs/scratch/(user
name)/petsc-3.6.2/src/mat/impls/aij/seq/aij.c
[0]PETSC ERROR: #5 MatCopy() line 3779 in /gpfs/scratch/(user
name)/petsc-3.6.2/src/mat/interface/matrix.c
48 MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
(gdb) n
49 MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
I tried a DMDACreate1D where I specify the dimension as n^3, where n is
the number of nodes on one side of the cube, but I still get the same
error. The error is generated at the MatCopy line as mentioned earlier.

This will not give you what you want since its got 1D connectivity.

Post by K. N. Ramachandran
In parallel you will not be able to reuse some old code to "compute the

Post by Barry Smith
Finite difference discretization of a 3D regular grid" so you might
consider just discarding that old code and using what PETSc provides for
computing the "Finite difference discretization of a 3D regular grid."

Ah ok, I had meant to say Finite difference discretization of the Poisson
equation.
If my parallel code could generate the sparse matrices formats for the
portion owned by each ranks, would there be a convenient way to generate
the PETSc matrices to be part of the DM? I had noticed that
CreateMatSeqAIJWithArrays changes the location of the Mat pointer passed in
for ComputeMatrix, so the actual matrix in DM remains as all zeros, so my
final solution comes out as zeros too.

This seems like the hardest way to do this. We have running examples, that
scale well, and produce exactly the matrix
you are using. In addition, they create the matrix in parallel, so the
whole thing is scalable.

In order to do it the way you want, you would need to match the
partitioning, which is hard.

Matt

Thanking You,

Post by K. N. Ramachandran
K.N.Ramachandran
Ph: 814-441-4279

[K. N. Ramachandran]

Hello Matt,

Post by Matthew Knepley
This seems like the hardest way to do this. We have running examples, that
scale well, and produce exactly the matrix
you are using. In addition, they create the matrix in parallel, so the
whole thing is scalable.
In order to do it the way you want, you would need to match the
partitioning, which is hard.

Yes, I understand. I am trying to keep the solver module separate from the
problem physics (like boundary conditions), i.e. my code has knowledge of
the physics and I generate the matrix in sparse format. I then hand it over
to some Solver module (say, PETSc) which can efficiently solve this.

Examples like snes/ex5.c, ksp/ex45.c, etc. use knowledge of the boundary
conditions to fill up the matrix and the RHS and this makes sense, if the
entire code is written using PETSc's native data structures. But from a
modularity point of view, if my problem needs a different partitioning
based on the physics, I feel like my code should generate the relevant
matrices on each rank and send it to PETSc.

I am trying to simulate the movement of charged particles in a domain, and
these particles become collimated along a particular direction. So one way
to partition might be 1D slabs perpendicular to the collimated direction,
so each rank can calculate the field and move the particles on each
timestep. Of course communication across the slabs would be inevitable. But
this way, the partitioning is based on the physics of the problem and the
solver module just efficiently solves the system without having to know the
boundary conditions, which makes sense to me overall.

Hope this description helps.


Thanking You,
K.N.Ramachandran
Ph: 814-441-4279

[Matthew Knepley]

Post by K. N. Ramachandran
Hello Matt,

Post by Matthew Knepley
This seems like the hardest way to do this. We have running examples,
that scale well, and produce exactly the matrix
you are using. In addition, they create the matrix in parallel, so the
whole thing is scalable.
In order to do it the way you want, you would need to match the
partitioning, which is hard.

Yes, I understand. I am trying to keep the solver module separate from the
problem physics (like boundary conditions), i.e. my code has knowledge of
the physics and I generate the matrix in sparse format. I then hand it over
to some Solver module (say, PETSc) which can efficiently solve this.
Examples like snes/ex5.c, ksp/ex45.c, etc. use knowledge of the boundary
conditions to fill up the matrix and the RHS and this makes sense, if the
entire code is written using PETSc's native data structures. But from a
modularity point of view, if my problem needs a different partitioning
based on the physics, I feel like my code should generate the relevant
matrices on each rank and send it to PETSc.
I am trying to simulate the movement of charged particles in a domain, and
these particles become collimated along a particular direction. So one way
to partition might be 1D slabs perpendicular to the collimated direction,
so each rank can calculate the field and move the particles on each
timestep. Of course communication across the slabs would be inevitable. But
this way, the partitioning is based on the physics of the problem and the
solver module just efficiently solves the system without having to know the
boundary conditions, which makes sense to me overall.
Hope this description helps.

The bad design decision is using a storage format designed for high speed
MatMult() to use for data interchange. You should be
using the API for value insertion rather than some storage format.

Matt

Post by K. N. Ramachandran
Thanking You,
K.N.Ramachandran
Ph: 814-441-4279

[K. N. Ramachandran]

Hello Matt,

Post by Matthew Knepley
The bad design decision is using a storage format designed for high speed
MatMult() to use for data interchange. You should be
using the API for value insertion rather than some storage format.

Ok, noted. Instead of using CreateMatSeqAIJWithArrays, I have also tried
setting up the Mat passed into the ComputeMatrix function, using
MatSetValue where the cols and values are referred to values in the CSR
format, but I still get the same error. It looks like the Mat is not
preallocated properly in the DM, even when I use DMDACreate3d.

Could you suggest what I can do to resolve this?

Thanks,
K.N.Ramachandran
Ph: 814-441-4279

[Matthew Knepley]

Post by K. N. Ramachandran
Hello Matt,

Post by Matthew Knepley
The bad design decision is using a storage format designed for high speed
MatMult() to use for data interchange. You should be
using the API for value insertion rather than some storage format.

Ok, noted. Instead of using CreateMatSeqAIJWithArrays, I have also tried
setting up the Mat passed into the ComputeMatrix function, using
MatSetValue where the cols and values are referred to values in the CSR
format, but I still get the same error. It looks like the Mat is not
preallocated properly in the DM, even when I use DMDACreate3d.

Another explanation is that the matrix is preallocated correctly, but you
are mapping the row and columns indices differently than the
DM. This what I was refering to when I said that getting the parallel
partition correct would be challenging.

Matt

Post by K. N. Ramachandran
Could you suggest what I can do to resolve this?
Thanks,
K.N.Ramachandran
Ph: 814-441-4279

[Jed Brown]

Post by K. N. Ramachandran
I am trying to simulate the movement of charged particles in a domain, and
these particles become collimated along a particular direction. So one way
to partition might be 1D slabs perpendicular to the collimated direction,
so each rank can calculate the field and move the particles on each
timestep. Of course communication across the slabs would be inevitable. But
this way, the partitioning is based on the physics of the problem and the
solver module just efficiently solves the system without having to know the
boundary conditions, which makes sense to me overall.

I don't understand what you mean by "without having to know the boundary
conditions". The boundary conditions are essential to solving
equations. Of course if you're assembling a matrix, they need to be
part of that matrix. (Usually the matrix will be singular "without"
boundary conditions.)

As for using a 1D partition, that may be very inefficient for the
solver, both due to limited scalability of the partition and more
interestingly if the directions of strong coupling do not lie entirely
within your planes and completely swamp the out-of-plane coupling
strength, then your solver can be expected to converge more slowly.

Note that DMDA can be instructed to produce a 1D partition if you like.

[Barry Smith]

Where is the matrix coming from? If it is from a 1, 2 or 3d structured grid then you can use the DMDA plus PCMG to manage almost everything including the generation of the interpolation. If you are not working with a structured grid then you can us PCGAMG which does algebraic multigrid using the matrix you provide.

Barry

Post by K. N. Ramachandran
Hello Matt,
Ah ok. I had tried specifying -pc_mg_levels earlier, but it asks for interpolation to be setup. I am currently working on using DMDA from examples 32 and 34, which sets up the interpolation, using DMDASetInterpolationType. So this is the recommended way to go forward right?
Since I supply the sparse matrix format and the rhs values to PETSc, does it make sense to do a DMDACreate1d instead of 3d? I did a sample implementation using DMDA 1d and when I call KSPSolve and fetch the solution, it comes out all zeros, with the solver having converged. I am trying to understand what I am doing wrong there.
Hello Matt,
Actually I felt the boundary conditions were having a role to play and set all the boundary conditions to Dirichlet. In this case, convergence was almost immediate with the Hypre preconditioner, taking 17 seconds with 3 iterations. The MG method took about the same time though.
So I reverted to the Dirichlet, Neumann mix of BCs and Hypre starts to diverge. Please find attached the output for the Hypre run using Dirichlet and Neumann for a 21^3 grid (rows), with a max of 7 nonzeros per row. Details of the options used before running are in the file. The solver used in all cases is bcgs.
Also attached is the MG output for 101^3 grid
1) Dirichlet and Neumann
2) Dirichlet only
where it seems to take about the same time.
Notice that you have no levels of MG here. You need to use -pc_mg_levels <n>
Matt
Thanking You,
Ramachandran K.N.
Ph: 814-441-4279

[Matthew Knepley]

Post by K. N. Ramachandran
Hello all,
I am a graduate student pursuing my Master's and I am trying to benchmark
a previous work by using PETSc for solving Poisson's Equation.
I am starting off with a serial code and I am trying to keep my code
modular, i.e. I generate the sparse matrix format and send it to PETSc or
any other solver. So I haven't built my code from the ground up using
PETSc's native data structures.
I am having trouble understanding the behavior of the solver and would
like your thoughts or inputs on what I can do better. I have both Dirichlet
and Neumann boundary conditions and my matrix size (number of rows) is
around a million but very sparse (~7 nonzeros per row), as can be expected
from a finite difference discretization of Poisson's equation.
I tried the methods outlined here
<http://scicomp.stackexchange.com/questions/513/why-is-my-iterative-linear-solver-not-converging?rq=1>
and here
<http://scicomp.stackexchange.com/questions/34/how-can-i-estimate-the-condition-number-of-a-large-sparse-matrix-using-petsc>.
Reverting to a 41^3 grid, I got the approximate condition number (using -ksp_monitor_singular_value
-ksp_type gmres -ksp_gmres_restart 1000 -pc_type none) as ~9072, which
seems pretty large. Higher matrix sizes give a larger condition number.
1) My best performing solver + preconditioner is bcgs+ilu(0) (on 1e6
grid) which solves in around 32 seconds, 196 iterations. How do I get a
fix for what the lower bound on the running time could be?

Unfortunately, that is really slow. You can look at this Tutorial (
http://www.mcs.anl.gov/petsc/documentation/tutorials/ParisTutorial.pdf)
slides 125-126 where
I solve elliptic systems using MG for systems with > 1M unknowns in that
amount of time.

You need to use Multigrid. Its the only solver that should be used for this
problem. Anything else is a waste of time.

Post by K. N. Ramachandran
2) Initially -pc_type hypre just Diverged and I was never able to use it.
Looking at this thread
<http://lists.mcs.anl.gov/pipermail/petsc-users/2013-October/019127.html>,
I had tried the options and it no longer diverges, but the residuals reduce
and hover around a constant value. How do I work with hypre to get a
useful preconditioner?

Something is wrong with your setup. Hypre is a fantastic solver for these
problems. Send the convergence output (-ksp_monitor_true_residual
-ksp_converged_reason -ksp_view).

Post by K. N. Ramachandran
Initially I solve Laplace's equation, so the mesh grid size has no effect
and even when I solve Poisson's equation, the spacing is carried over to
the RHS, so I am pretty sure the spacing is not affecting the condition
number calculation.

No, the Laplacian has a condition number that grows like h^{-2}. See
earlier slides in that tutorial.

Thanks,

Matt

Post by K. N. Ramachandran
Hope this helps. Please let me know if you might need more information.
Thanking You,
K.N.Ramachandran
Ph: 814-441-4279