Re: [petsc-users] Problems about Block Jocabi and Matrix-Free Method in SNES

2018-10-10 Thread Matthew Knepley
On Wed, Oct 10, 2018 at 8:50 AM Yingjie Wu  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,);
> 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:   Vect;
> 289:   PetscInt   m;
>
> 292:   DMGetGlobalVector(packer,);
> 293:   VecGetLocalSize(t,);
> 294:   DMRestoreGlobalVector(packer,);
> 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*));
>
> 131:   PetscOptionsGetBool(NULL,NULL,"-snes_mf",_free,NULL);
> 132:   PetscOptionsGetBool(NULL,NULL,"-fdcoloring",,NULL);
> 133:   PetscOptionsGetBool(NULL,NULL,"-fdcoloring_ds",_ds,NULL);
> 134:
>  PetscOptionsGetBool(NULL,NULL,"-fdcoloring_local",_coloring,NULL);
> 135:   if (!matrix_free) {
> 136: DMSetMatType(user.da,MATAIJ);
> 137: DMCreateMatrix(user.da,);
> 138: if (coloring) {
> 139:   ISColoring iscoloring;
> 140:   if (!local_coloring) {
> 141: DMCreateColoring(user.da,IS_COLORING_GLOBAL,);
> 142: MatFDColoringCreate(J,iscoloring,);
> 143: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode
> (*)(void))FormFunction,);
> 144:   } else {
> 145: DMCreateColoring(user.da,IS_COLORING_LOCAL,);
> 146: MatFDColoringCreate(J,iscoloring,);
> 147: MatFDColoringUseDM(J,matfdcoloring);
> 148: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode
> (*)(void))FormFunctionLocal,);
> 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();
> 157: } else {
> 158:   SNESSetJacobian(snes,J,J,FormJacobian,);
> 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  于2018年10月9日周二 下午11:34写道:
>
>> On Mon, Oct 8, 2018 at 11:33 PM Yingjie Wu  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 1 -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 1 -ksp_type fgmres -pc_type 

Re: [petsc-users] Problems about Block Jocabi and Matrix-Free Method in SNES

2018-10-10 Thread Yingjie Wu
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,);
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:   Vect;
289:   PetscInt   m;

292:   DMGetGlobalVector(packer,);
293:   VecGetLocalSize(t,);
294:   DMRestoreGlobalVector(packer,);
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*));

131:   PetscOptionsGetBool(NULL,NULL,"-snes_mf",_free,NULL);
132:   PetscOptionsGetBool(NULL,NULL,"-fdcoloring",,NULL);
133:   PetscOptionsGetBool(NULL,NULL,"-fdcoloring_ds",_ds,NULL);
134:
 PetscOptionsGetBool(NULL,NULL,"-fdcoloring_local",_coloring,NULL);
135:   if (!matrix_free) {
136: DMSetMatType(user.da,MATAIJ);
137: DMCreateMatrix(user.da,);
138: if (coloring) {
139:   ISColoring iscoloring;
140:   if (!local_coloring) {
141: DMCreateColoring(user.da,IS_COLORING_GLOBAL,);
142: MatFDColoringCreate(J,iscoloring,);
143: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode
(*)(void))FormFunction,);
144:   } else {
145: DMCreateColoring(user.da,IS_COLORING_LOCAL,);
146: MatFDColoringCreate(J,iscoloring,);
147: MatFDColoringUseDM(J,matfdcoloring);
148: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode
(*)(void))FormFunctionLocal,);
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();
157: } else {
158:   SNESSetJacobian(snes,J,J,FormJacobian,);
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  于2018年10月9日周二 下午11:34写道:

> On Mon, Oct 8, 2018 at 11:33 PM Yingjie Wu  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 1 -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 1 -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 

Re: [petsc-users] Problems about Block Jocabi and Matrix-Free Method in SNES

2018-10-09 Thread Matthew Knepley
On Mon, Oct 8, 2018 at 11:33 PM Yingjie Wu  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 1 -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 1 -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 

[petsc-users] Problems about Block Jocabi and Matrix-Free Method in SNES

2018-10-08 Thread Yingjie Wu
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 1 -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
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
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