On Wed, Oct 10, 2018 at 8:50 AM Yingjie Wu <yjw...@gmail.com> wrote: > Thank you very much, Matt. > I followed your advice and learned some MATMFFD related functions. Since > I really wanted to use the Matrix-Free method, I learned the few examples > of this method. I would like to enhance my understanding of the Matrix-Free > method in PETSc according to my questions about these examples and hope to > get your answer. > 1. In the example src/snes/example/tutorials/ex22.c > > ...... > 107: packer->ops->creatematrix = DMCreateMatrix_MF; > > 110: SNESCreate(PETSC_COMM_WORLD,&snes); > 111: SNESSetDM(snes,packer); > 112: SNESSetFunction(snes,NULL,ComputeFunction,NULL); > 113: SNESSetJacobian(snes,NULL, NULL,ComputeJacobian_MF,NULL); > > 115: SNESSetFromOptions(snes); > ...... > 285: PetscErrorCode DMCreateMatrix_MF(DM packer,Mat *A) > 286: { > 288: Vec t; > 289: PetscInt m; > > 292: DMGetGlobalVector(packer,&t); > 293: VecGetLocalSize(t,&m); > 294: DMRestoreGlobalVector(packer,&t); > 295: > MatCreateMFFD(PETSC_COMM_WORLD,m,m,PETSC_DETERMINE,PETSC_DETERMINE,A); > 296: MatSetUp(*A); > 297: return(0); > 298: } > > 300: PetscErrorCode ComputeJacobian_MF(SNES snes,Vec x,Mat A,Mat B,void > *ctx) > 301: { > > 305: MatMFFDSetFunction(A,(PetscErrorCode > (*)(void*,Vec,Vec))SNESComputeFunction,snes); > 306: MatMFFDSetBase(A,x,NULL); > 307: return(0); > 308: } > > Why are the solution vector, Jacobian matrix, and precondition matrix set > to NULL in SNESSetJacobian? Is it because they are automatically generated > in SNESSetDM? >
Yes. This example show how to do it by hand. If you give NULL for your FormJacobian function, we do it automatically. If you give _mf_operator, then we do it, but pass through the matrices to your FormJacobian() function. > There is no precondition matrix B in ComputeJacobian_MF (Jacobian evaluate > program). If I want to add my own precondition matrix, then the Assembly > matrix B in this function, is that right? > What is the function of MatMFFDSetBase? I read the documentation but > didn't understand it. What would happen if I removed it? > SetBase() provides the initial vector for differencing. This is necessary. > 2. In the example /src/snes/examples/tutorials/ex14.c > > 114: SNESSetFunction(snes,r,FormFunction,(void*)&user); > > 131: PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL); > 132: PetscOptionsGetBool(NULL,NULL,"-fdcoloring",&coloring,NULL); > 133: PetscOptionsGetBool(NULL,NULL,"-fdcoloring_ds",&coloring_ds,NULL); > 134: > PetscOptionsGetBool(NULL,NULL,"-fdcoloring_local",&local_coloring,NULL); > 135: if (!matrix_free) { > 136: DMSetMatType(user.da,MATAIJ); > 137: DMCreateMatrix(user.da,&J); > 138: if (coloring) { > 139: ISColoring iscoloring; > 140: if (!local_coloring) { > 141: DMCreateColoring(user.da,IS_COLORING_GLOBAL,&iscoloring); > 142: MatFDColoringCreate(J,iscoloring,&matfdcoloring); > 143: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode > (*)(void))FormFunction,&user); > 144: } else { > 145: DMCreateColoring(user.da,IS_COLORING_LOCAL,&iscoloring); > 146: MatFDColoringCreate(J,iscoloring,&matfdcoloring); > 147: MatFDColoringUseDM(J,matfdcoloring); > 148: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode > (*)(void))FormFunctionLocal,&user); > 149: } > 150: if (coloring_ds) { > 151: MatFDColoringSetType(matfdcoloring,MATMFFD_DS); > 152: } > 153: MatFDColoringSetFromOptions(matfdcoloring); > 154: MatFDColoringSetUp(J,iscoloring,matfdcoloring); > 155: > SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring); > 156: ISColoringDestroy(&iscoloring); > 157: } else { > 158: SNESSetJacobian(snes,J,J,FormJacobian,&user); > 159: } > 160: } > > If the command "- snes_mf" is used, the function SNESSetJacobian is not > called, and the Jacobian matrix type is initialized as MATMFFD. The > function in SNESSetFunction is used to approximate the matrix vector > product. Do I understand this correctly? > Yes . Thanks, Matt > The description of my question may be complicated, but as a beginner, I > really want to learn something about matrix-free method from the examples, > and I hope to get your reply. > > Thanks, > Yingjie > > Matthew Knepley <knep...@gmail.com> 于2018年10月9日周二 下午11:34写道: > >> On Mon, Oct 8, 2018 at 11:33 PM Yingjie Wu <yjw...@gmail.com> wrote: >> >>> Dear Petsc developer: >>> Hi, >>> >>> I've been studying Petsc recently about Precontioner and Metrix-Free, >>> and I have some questions that puzzle me. >>> >>> 1. I want to test block Jacobi preconditioner, so I chose >>> /snes/example/tutorial/ex3.c as an example. According to the reference in >>> the example, the input parameters are: >>> mpiexec -n 8./ex3 -nox -n 10000 -ksp_type fgmres -pc_type bjacobi >>> -pc_bjacobi_blocks 4 -sub_ksp_type gmres -sub_ksp_max_it 3 -post_setsubksp >>> -sub_ksp_rtol 1.e-16 >>> >> >> You do not care about recursive blocks, so just use >> >> $MPIEXEC -n 8 ./ex3 -nox -n 10000 -ksp_type fgmres -pc_type bjacobi >> -pc_bjacobi_blocks 4 -sub_ksp_type gmres -sub_ksp_max_it 3 -snes_view >> -ksp_view >> >> and I get the attached output. >> >> >>> I want to export each block of KSP and PC information : >>> -snes_view -ksp_view >>> However, the procedure is wrong and the wrong information is as >>> follows: >>> >>> [0]PETSC ERROR: --------------------- Error Message >>> -------------------------------------------------------------- >>> [0]PETSC ERROR: Arguments must have same communicators >>> [0]PETSC ERROR: Different communicators in the two objects: Argument # 1 >>> and 2 flag 3 >>> [0]PETSC ERROR: [1]PETSC ERROR: --------------------- Error Message >>> -------------------------------------------------------------- >>> [2]PETSC ERROR: --------------------- Error Message >>> -------------------------------------------------------------- >>> [2]PETSC ERROR: [3]PETSC ERROR: --------------------- Error Message >>> -------------------------------------------------------------- >>> [3]PETSC ERROR: Arguments must have same communicators >>> [3]PETSC ERROR: [6]PETSC ERROR: --------------------- Error Message >>> -------------------------------------------------------------- >>> [6]PETSC ERROR: Arguments must have same communicators >>> [6]PETSC ERROR: Different communicators in the two objects: Argument # 1 >>> and 2 flag 3 >>> [6]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html >>> for trouble shooting. >>> [7]PETSC ERROR: --------------------- Error Message >>> -------------------------------------------------------------- >>> [7]PETSC ERROR: Arguments must have same communicators >>> [7]PETSC ERROR: Different communicators in the two objects: Argument # 1 >>> and 2 flag 3 >>> [7]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html >>> for trouble shooting. >>> [7]PETSC ERROR: Petsc Release Version 3.10.1, Sep, 26, 2018 >>> See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble >>> shooting. >>> [0]PETSC ERROR: Petsc Release Version 3.10.1, Sep, 26, 2018 >>> [0]PETSC ERROR: ./ex3 on a arch-linux2-c-debug named yjwu-XPS-8910 by >>> yjwu Mon Oct 8 22:35:34 2018 >>> [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ >>> --with-fc=gfortran --download-mpich --download-fblaslapack >>> [0]PETSC ERROR: #1 KSPView() line 118 in >>> /home/yjwu/petsc-3.10.1/src/ksp/ksp/interface/itcreate.c >>> [0]PETSC ERROR: #2 PCView_BJacobi() line 232 in >>> /home/yjwu/petsc-3.10.1/src/ksp/pc/impls/bjacobi/bjacobi.c >>> [0]PETSC ERROR: #3 PCView() line 1651 in >>> /home/yjwu/petsc-3.10.1/src/ksp/pc/interface/precon.c >>> [1]PETSC ERROR: Arguments must have same communicators >>> [1]PETSC ERROR: Different communicators in the two objects: Argument # 1 >>> and 2 flag 3 >>> [1]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html >>> for trouble shooting. >>> [1]PETSC ERROR: Petsc Release Version 3.10.1, Sep, 26, 2018 >>> [1]PETSC ERROR: ./ex3 on a arch-linux2-c-debug named yjwu-XPS-8910 by >>> yjwu Mon Oct 8 22:35:34 2018 >>> [1]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ >>> --with-fc=gfortran --download-mpich --download-fblaslapack >>> [1]PETSC ERROR: #1 KSPView() line 118 in >>> /home/yjwu/petsc-3.10.1/src/ksp/ksp/interface/itcreate.c >>> Arguments must have same communicators >>> [2]PETSC ERROR: Different communicators in the two objects: Argument # 1 >>> and 2 flag 3 >>> [2]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html >>> for trouble shooting. >>> >>> If I want to get the information of KSP and PC in each block, what >>> should I do? >>> >>> 2. There is a requirement in my program to use Matrix-Free method to >>> approximate Jacobian matrix by finite difference of residual function in >>> solving nonlinear equations. But I'll also provide an analytic( incomplete, >>> some terms are missing or approximate) for preconditioning. Because my >>> problem is about a large set of equations composed of several physical >>> fields, I want to use block Jacobian precondition for each subfield(block), >>> and ILU sub_pc for each subfield. After reading the Users'Guide, I found >>> that using - snex_mf_operator can do the above, so I added: >>> - snes_mf_operator >>> >> >> You need to be careful what matrix you are adding values to in your >> FormJacobian() routine. The primal matrix J is of type MFFD (finite >> difference) >> and thus cannot accept values. You put your approximate values in the >> preconditioner M. >> >> Thanks, >> >> Matt >> >> >>> after the example above. However, the procedure is wrong and does not >>> seem to support this. >>> The wrong information is: >>> >>> [1]PETSC ERROR: --------------------- Error Message >>> -------------------------------------------------------------- >>> [2]PETSC ERROR: --------------------- Error Message >>> -------------------------------------------------------------- >>> [2]PETSC ERROR: No support for this operation for this object type >>> [2]PETSC ERROR: Mat type mffd >>> [2]PETSC ERROR: [3]PETSC ERROR: [4]PETSC ERROR: --------------------- >>> Error Message -------------------------------------------------------------- >>> [4]PETSC ERROR: No support for this operation for this object type >>> [4]PETSC ERROR: Mat type mffd >>> [4]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html >>> for trouble shooting. >>> [4]PETSC ERROR: Petsc Release Version 3.10.1, Sep, 26, 2018 >>> [6]PETSC ERROR: --------------------- Error Message >>> -------------------------------------------------------------- >>> [6]PETSC ERROR: No support for this operation for this object type >>> [6]PETSC ERROR: Mat type mffd >>> [6]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html >>> for trouble shooting. >>> [6]PETSC ERROR: Petsc Release Version 3.10.1, Sep, 26, 2018 >>> [6]PETSC ERROR: ./ex3 on a arch-linux2-c-debug named yjwu-XPS-8910 by >>> yjwu Mon Oct 8 23:01:57 2018 >>> [6]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ >>> --with-fc=gfortran --download-mpich --download-fblaslapack >>> [0]PETSC ERROR: --------------------- Error Message >>> -------------------------------------------------------------- >>> [0]PETSC ERROR: No support for this operation for this object type >>> [0]PETSC ERROR: Mat type mffd >>> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html >>> for trouble shooting. >>> [0]PETSC ERROR: Petsc Release Version 3.10.1, Sep, 26, 2018 >>> [0]PETSC ERROR: ./ex3 on a arch-linux2-c-debug named yjwu-XPS-8910 by >>> yjwu Mon Oct 8 23:01:57 2018 >>> [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ >>> --with-fc=gfortran --download-mpich --download-fblaslapack >>> [0]PETSC ERROR: #1 MatSetValues() line 1330 in >>> /home/yjwu/petsc-3.10.1/src/mat/interface/matrix.c >>> [1]PETSC ERROR: No support for this operation for this object type >>> [1]PETSC ERROR: Mat type mffd >>> [1]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html >>> for trouble shooting. >>> >>> If I want to implement matrix-free and block jacobi precondition in my >>> program, what do I need to pay attention or add ? >>> >>> The second question is very important for my research and I hope to get >>> your answer. >>> >>> Thanks, >>> Yingjie >>> >> >> >> -- >> 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 >> >> https://www.cse.buffalo.edu/~knepley/ >> <http://www.cse.buffalo.edu/~knepley/> >> > -- 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 https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>