[petsc-users] SNESComputeJacobianDefaultColor troubles

Discussion:

Francesco Magaletti

2015-11-02 15:58:51 UTC

Hi everyone,

I’m trying to solve a system of PDE’s with a full implicit time integration and a DMDA context to manage the cartesian structured grid.
I don’t know the actual jacobian matrix (actually it is too cumbersome to be easily evaluated analytically) so I used SNESComputeJacobianDefaultColor to approximate it:

DMSetMatType(da,MATAIJ);
DMCreateMatrix(da,&J);
TSGetSNES(ts, &snes);
SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor,0);

All it works perfectly but it happens that, when I increase the number of grid points, I go out of memory since the matrix becomes too big to be stored despite my estimate of memory consumption which is much lower than the actual one. Indeed a more careful analysis (with -mat_view) shows that there are a lot of zeros in the nonzero structure of the jacobian. I read in the manual that DMCreateMatrix preallocates the matrix by somewhat using the stencil width of the DMDA context, but in my case I don’t really need all those elements for every equation.
I then tried with a poor preallocation, hoping petsC could manage it with some malloc on the fly:

MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,size,size,3,0,0,0,&J);
MatSetOption(J,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
TSGetSNES(ts, &snes);
SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor,0);

but now it gives me runtime errors:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: Matrix must be assembled by calls to MatAssemblyBegin/End();
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015
[0]PETSC ERROR: ./bubble_petsc on a arch-linux2-c-debug named ... Mon Nov 2 16:24:12 2015
[0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich
[0]PETSC ERROR: #1 MatFDColoringCreate() line 448 in /home/magaletto/Scaricati/petsc-3.6.0/src/mat/matfd/fdmatrix.c
[0]PETSC ERROR: #2 SNESComputeJacobianDefaultColor() line 67 in /home/magaletto/Scaricati/petsc-3.6.0/src/snes/interface/snesj2.c
[0]PETSC ERROR: #3 SNESComputeJacobian() line 2223 in /home/magaletto/Scaricati/petsc-3.6.0/src/snes/interface/snes.c
[0]PETSC ERROR: #4 SNESSolve_NEWTONLS() line 231 in /home/magaletto/Scaricati/petsc-3.6.0/src/snes/impls/ls/ls.c
[0]PETSC ERROR: #5 SNESSolve() line 3894 in /home/magaletto/Scaricati/petsc-3.6.0/src/snes/interface/snes.c
[0]PETSC ERROR: #6 TSStep_Theta() line 197 in /home/magaletto/Scaricati/petsc-3.6.0/src/ts/impls/implicit/theta/theta.c
[0]PETSC ERROR: #7 TSStep() line 3098 in /home/magaletto/Scaricati/petsc-3.6.0/src/ts/interface/ts.c
[0]PETSC ERROR: #8 TSSolve() line 3282 in /home/magaletto/Scaricati/petsc-3.6.0/src/ts/interface/ts.c

If I use instead the same preallocation but with the FD calculation without coloring it runs without problems (it is only extremely slow) and it uses the correct number of nonzero values:

MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,size,size,3,0,0,0,&J);
MatSetOption(J,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
TSGetSNES(ts, &snes);
SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault,&appctx);

What do you suggest to solve this problem?
I also tried to set up a first step with fd without coloring in order to allocate the correct amount of memory and nonzero structure and successively to switch to the colored version:

TSSetIFunction(ts,NULL,FormIFunction,&appctx);

MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,size,size,3,0,0,0,&J);
MatSetOption(J,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
TSGetSNES(ts, &snes);
SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault,&appctx);
TSStep(ts);

ISColoring iscoloring;
MatColoring coloring;
MatFDColoring matfdcoloring;
MatColoringCreate(J,&coloring);
MatColoringSetFromOptions(coloring);
MatColoringApply(coloring,&iscoloring);
MatFDColoringCreate(J,iscoloring,&matfdcoloring);
MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormIFunction,&appctx);
MatFDColoringSetFromOptions(matfdcoloring);
MatFDColoringSetUp(J,iscoloring,matfdcoloring);
SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);

TSStep(ts);

but again it gives me runtime errors:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Invalid argument
[0]PETSC ERROR: Wrong type of object: Parameter # 1
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015
[0]PETSC ERROR: ./bubble_petsc on a arch-linux2-c-debug named ... Mon Nov 2 16:37:24 2015
[0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich
[0]PETSC ERROR: #1 TSGetDM() line 4093 in /home/magaletto/Scaricati/petsc-3.6.0/src/ts/interface/ts.c
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Corrupt argument: http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
[0]PETSC ERROR: Invalid Pointer to Object: Parameter # 1
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015
[0]PETSC ERROR: ./bubble_petsc on a arch-linux2-c-debug named ... Mon Nov 2 16:37:24 2015
[0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich
[0]PETSC ERROR: #2 DMGetLocalVector() line 44 in /home/magaletto/Scaricati/petsc-3.6.0/src/dm/interface/dmget.c

=====================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= EXIT CODE: 11
= CLEANING UP REMAINING PROCESSES
= YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
=====================================================================================
APPLICATION TERMINATED WITH THE EXIT STRING: Segmentation fault (signal 11)

where is my mistake?

Your sincerely,
Francesco

Barry Smith

2015-11-02 16:54:00 UTC

Permalink

Don't do any of that stuff you tried below. Just use DMDASetFillBlocks() before DMCreateMatrix()

Barry

Post by Francesco Magaletti
Hi everyone,
I’m trying to solve a system of PDE’s with a full implicit time integration and a DMDA context to manage the cartesian structured grid.
DMSetMatType(da,MATAIJ);
DMCreateMatrix(da,&J);
TSGetSNES(ts, &snes);
SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor,0);
All it works perfectly but it happens that, when I increase the number of grid points, I go out of memory since the matrix becomes too big to be stored despite my estimate of memory consumption which is much lower than the actual one. Indeed a more careful analysis (with -mat_view) shows that there are a lot of zeros in the nonzero structure of the jacobian. I read in the manual that DMCreateMatrix preallocates the matrix by somewhat using the stencil width of the DMDA context, but in my case I don’t really need all those elements for every equation.
MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,size,size,3,0,0,0,&J);
MatSetOption(J,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
TSGetSNES(ts, &snes);
SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor,0);
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: Matrix must be assembled by calls to MatAssemblyBegin/End();
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015
[0]PETSC ERROR: ./bubble_petsc on a arch-linux2-c-debug named ... Mon Nov 2 16:24:12 2015
[0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich
[0]PETSC ERROR: #1 MatFDColoringCreate() line 448 in /home/magaletto/Scaricati/petsc-3.6.0/src/mat/matfd/fdmatrix.c
[0]PETSC ERROR: #2 SNESComputeJacobianDefaultColor() line 67 in /home/magaletto/Scaricati/petsc-3.6.0/src/snes/interface/snesj2.c
[0]PETSC ERROR: #3 SNESComputeJacobian() line 2223 in /home/magaletto/Scaricati/petsc-3.6.0/src/snes/interface/snes.c
[0]PETSC ERROR: #4 SNESSolve_NEWTONLS() line 231 in /home/magaletto/Scaricati/petsc-3.6.0/src/snes/impls/ls/ls.c
[0]PETSC ERROR: #5 SNESSolve() line 3894 in /home/magaletto/Scaricati/petsc-3.6.0/src/snes/interface/snes.c
[0]PETSC ERROR: #6 TSStep_Theta() line 197 in /home/magaletto/Scaricati/petsc-3.6.0/src/ts/impls/implicit/theta/theta.c
[0]PETSC ERROR: #7 TSStep() line 3098 in /home/magaletto/Scaricati/petsc-3.6.0/src/ts/interface/ts.c
[0]PETSC ERROR: #8 TSSolve() line 3282 in /home/magaletto/Scaricati/petsc-3.6.0/src/ts/interface/ts.c
MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,size,size,3,0,0,0,&J);
MatSetOption(J,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
TSGetSNES(ts, &snes);
SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault,&appctx);
What do you suggest to solve this problem?
TSSetIFunction(ts,NULL,FormIFunction,&appctx);
MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,size,size,3,0,0,0,&J);
MatSetOption(J,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
TSGetSNES(ts, &snes);
SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault,&appctx);
TSStep(ts);
ISColoring iscoloring;
MatColoring coloring;
MatFDColoring matfdcoloring;
MatColoringCreate(J,&coloring);
MatColoringSetFromOptions(coloring);
MatColoringApply(coloring,&iscoloring);
MatFDColoringCreate(J,iscoloring,&matfdcoloring);
MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormIFunction,&appctx);
MatFDColoringSetFromOptions(matfdcoloring);
MatFDColoringSetUp(J,iscoloring,matfdcoloring);
SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);
TSStep(ts);
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Invalid argument
[0]PETSC ERROR: Wrong type of object: Parameter # 1
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015
[0]PETSC ERROR: ./bubble_petsc on a arch-linux2-c-debug named ... Mon Nov 2 16:37:24 2015
[0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich
[0]PETSC ERROR: #1 TSGetDM() line 4093 in /home/magaletto/Scaricati/petsc-3.6.0/src/ts/interface/ts.c
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Corrupt argument: http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
[0]PETSC ERROR: Invalid Pointer to Object: Parameter # 1
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015
[0]PETSC ERROR: ./bubble_petsc on a arch-linux2-c-debug named ... Mon Nov 2 16:37:24 2015
[0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich
[0]PETSC ERROR: #2 DMGetLocalVector() line 44 in /home/magaletto/Scaricati/petsc-3.6.0/src/dm/interface/dmget.c
=====================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= EXIT CODE: 11
= CLEANING UP REMAINING PROCESSES
= YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
=====================================================================================
APPLICATION TERMINATED WITH THE EXIT STRING: Segmentation fault (signal 11)
where is my mistake?
Your sincerely,
Francesco

Francesco Magaletti

2015-11-02 22:06:11 UTC

Permalink

Thanks Barry!!
I tried and it worked nice, but it will be improved if I could manage the off-diagonal elements block by block.
I have a stencil with 5 points (in 1D) and the coupling with the i+1 is different with respect to the i+2 elements. I found on PetsC site a reference to DMDASetGetMatrix function, but there are no examples. Can you give me an example with a simple matrix? Or there is a simpler way to treat the different off-diagonal blocks separately?

Thank you again
Francesco

Post by Barry Smith
Don't do any of that stuff you tried below. Just use DMDASetFillBlocks() before DMCreateMatrix()
Barry

Barry Smith

2015-11-02 22:26:46 UTC

Permalink

Post by Francesco Magaletti
Thanks Barry!!
I tried and it worked nice, but it will be improved if I could manage the off-diagonal elements block by block.
I have a stencil with 5 points (in 1D) and the coupling with the i+1 is different with respect to the i+2 elements. I found on PetsC site a reference to DMDASetGetMatrix function, but there are no examples. Can you give me an example with a simple matrix? Or there is a simpler way to treat the different off-diagonal blocks separately?

You will need to dig directly into the PETSc code that generates the matrix preallocation and nonzero structure for the DMDA and make a copy of the routine and modify it as needed. See for example DMCreateMatrix_DA_2d_MPIAIJ_Fill in src/dm/impls/da/fdda.c

Post by Francesco Magaletti
Thank you again
Francesco

Post by Barry Smith
Don't do any of that stuff you tried below. Just use DMDASetFillBlocks() before DMCreateMatrix()
Barry