[petsc-users] Problems about Block Jocabi and Matrix-Free Method in SNES
Yingjie Wu
yjwu16 at gmail.com
Wed Oct 10 07:50:20 CDT 2018
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?
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?
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?
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/>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20181010/163a36ab/attachment-0001.html>
More information about the petsc-users
mailing list