[petsc-users] Problems about Block Jocabi and Matrix-Free Method in SNES
Matthew Knepley
knepley at gmail.com
Wed Oct 10 08:09:29 CDT 2018
On Wed, Oct 10, 2018 at 8:50 AM Yingjie Wu <yjwu16 at 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 <knepley at gmail.com> 于2018年10月9日周二 下午11:34写道：
>
>> On Mon, Oct 8, 2018 at 11:33 PM Yingjie Wu <yjwu16 at 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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20181010/d7e0fa60/attachment-0001.html>
More information about the petsc-users
mailing list