> On Nov 2, 2015, at 4:06 PM, Francesco Magaletti > <[email protected]> wrote: > > 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 > > Thank you again > Francesco > >> Il giorno 02/nov/2015, alle ore 17:54, Barry Smith <[email protected]> ha >> scritto: >> >> >> Don't do any of that stuff you tried below. Just use DMDASetFillBlocks() >> before DMCreateMatrix() >> >> Barry >> >> >>> On Nov 2, 2015, at 9:58 AM, Francesco Magaletti >>> <[email protected]> wrote: >>> >>> 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 >> >
