[petsc-users] Problems about Assemble DMComposite Precondition Matrix

[Jed Brown]

Yingjie Wu via petsc-users <petsc-users at mcs.anl.gov> writes:

> Thank you very much for your reply.
> My equation is a neutron diffusion equation with eigenvalues, which is why
> I use DMConposite because there is a single non-physical field variable,
> eigenvalue. I am not very familiar with FieldSplit. I can understand it
> first.
> It seems not the problem of DmConposite because the reason of error
> reporting is that the diagonal elements of the precondition matrix have
> zero terms.

Which entries are zero?  You can view the matrix to find out if it has
the structure you intend.  Make sure the index sets is[] are what you
intend.

> My requirement is to divide precondition matrix to sub-matrices, because I
> have neutrons of multiple energy groups (each is a physical field). I want
> to assign only diagonal block matrices, preferably using
> MatSetValuesStencil (which can simplify the assignment of five diagonal
> matrices).

Splitting all of those apart with DMComposite is clumsy and may deliver
poor memory performance.  Consider DMDASetBlockFills if you want to
allocate a matrix that has sparsity within blocks.

>
> Thanks,
> Yingjie
>
> Mark Adams <mfadams at lbl.gov> 于2018年11月5日周一 下午11:22写道：
>
>> DMComposite is not very mature, the last time I checked and I don't of
>> anyone having worked on it recently, and it is probably not what you want
>> anyway. FieldSplit is most likely what you want.
>>
>> What are your equations and discretization? eg, Stokes with cell centered
>> pressure? There are probably examples that are close to what you want and
>> it would not be hard to move your code over.
>>
>> Mark
>>
>> On Mon, Nov 5, 2018 at 10:00 AM Yingjie Wu via petsc-users <
>> petsc-users at mcs.anl.gov> wrote:
>>
>>> Dear Petsc developer:
>>> Hi,
>>> I have recently studied the preconditioner of the program, and some
>>> problems have arisen. Please help me to solve them.
>>> At present, I have written a program to solve the system of non-linear
>>> equations. The Matrix Free method has been used to calculate results. But I
>>> want to add a preprocessing matrix to it.
>>> I used the DMComposite object, which stores two sub-DM objects and a
>>> single value (two physical field variables and one variable). I want to use
>>> MatGetLocalSubMatrix to assign each physical field sub precondition matrix.
>>> At the same time, because my DM object is two-dimensional, I use
>>> MatSetValuesStencil() to fill the sub matrix.
>>> At present, I just want to fill in a unit matrix (for global vectors) as
>>> the precondition matrix of Matrix Free Method (just to test whether it can
>>> be used, the unit matrix has no preprocessing effect). But the procedure
>>> was wrong.
>>>
>>>            yjwu at yjwu-XPS-8910:~/petsc-3.10.1/src/snes/examples/tutorials$
>>> mpiexec -n 1 ./ex216 -f wu-readtwogroups -snes_mf_operator -snes_view
>>> -snes_converged_reason -snes_monitor -ksp_converged_reason
>>> -ksp_monitor_true_residual
>>>
>>>   0 SNES Function norm 8.235090086536e-02
>>> iter = 0, SNES Function norm 0.0823509
>>> iter = 0, Keff ======= 1.
>>>   Linear solve did not converge due to DIVERGED_PCSETUP_FAILED iterations
>>> 0
>>>                  PCSETUP_FAILED due to FACTOR_NUMERIC_ZEROPIVOT
>>> Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
>>> SNES Object: 1 MPI processes
>>>   type: newtonls
>>>   maximum iterations=50, maximum function evaluations=10000
>>>   tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
>>>   total number of linear solver iterations=0
>>>   total number of function evaluations=1
>>>   norm schedule ALWAYS
>>>   SNESLineSearch Object: 1 MPI processes
>>>     type: bt
>>>       interpolation: cubic
>>>       alpha=1.000000e-04
>>>     maxstep=1.000000e+08, minlambda=1.000000e-12
>>>     tolerances: relative=1.000000e-08, absolute=1.000000e-15,
>>> lambda=1.000000e-08
>>>     maximum iterations=40
>>>   KSP Object: 1 MPI processes
>>>     type: gmres
>>>       restart=30, using Classical (unmodified) Gram-Schmidt
>>> Orthogonalization with no iterative refinement
>>>       happy breakdown tolerance 1e-30
>>>     maximum iterations=10000, initial guess is zero
>>>     tolerances:  relative=1e-05, absolute=1e-50, divergence=10000.
>>>     left preconditioning
>>>     using PRECONDITIONED norm type for convergence test
>>>   PC Object: 1 MPI processes
>>>     type: ilu
>>>       out-of-place factorization
>>>       0 levels of fill
>>>       tolerance for zero pivot 2.22045e-14
>>>       matrix ordering: natural
>>>       factor fill ratio given 1., needed 1.
>>>         Factored matrix follows:
>>>           Mat Object: 1 MPI processes
>>>             type: seqaij
>>>             rows=961, cols=961
>>>             package used to perform factorization: petsc
>>>             total: nonzeros=4625, allocated nonzeros=4625
>>>             total number of mallocs used during MatSetValues calls =0
>>>               not using I-node routines
>>>     linear system matrix followed by preconditioner matrix:
>>>     Mat Object: 1 MPI processes
>>>       type: mffd
>>>       rows=961, cols=961
>>>         Matrix-free approximation:
>>>           err=1.49012e-08 (relative error in function evaluation)
>>>           The compute h routine has not yet been set
>>>     Mat Object: 1 MPI processes
>>>       type: seqaij
>>>       rows=961, cols=961
>>>       total: nonzeros=4625, allocated nonzeros=4625
>>>       total number of mallocs used during MatSetValues calls =0
>>>         not using I-node routines
>>>
>>> It seems that there are elements on the diagonal line that are not
>>> filled, but I don't understand what went wrong. Part of the code of my
>>> program is as follows. I added all the code to the appendix. I have tested
>>> that for the entire Jacobian matrix assignment unit matrix (no submatrix,
>>> direct operation of the global preprocessing matrix), the program can run
>>> normally. However, since the problems developed later may be more complex,
>>> involving multiple physical fields, it would be easier to implement
>>> submatrix.
>>>
>>> ......
>>>
>>> /* set DMComposite */
>>>
>>> ierr =
>>> DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,20,24,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&user.da1);CHKERRQ(ierr);
>>>   ierr = DMSetFromOptions(user.da1);CHKERRQ(ierr);
>>>   ierr = DMSetUp(user.da1);CHKERRQ(ierr);
>>>   ierr = DMCompositeAddDM(user.packer,user.da1);CHKERRQ(ierr);
>>>   ierr =
>>> DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,20,24,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&user.da2);CHKERRQ(ierr);
>>>   ierr = DMSetFromOptions(user.da2);CHKERRQ(ierr);
>>>   ierr = DMSetUp(user.da2);CHKERRQ(ierr);
>>>   ierr = DMCompositeAddDM(user.packer,user.da2);CHKERRQ(ierr);
>>>   ierr = DMRedundantCreate(PETSC_COMM_WORLD,0,1,&user.red1);CHKERRQ(ierr);
>>>   ierr = DMCompositeAddDM(user.packer,user.red1);CHKERRQ(ierr);
>>> ......
>>> /* in FormJacobian(SNES snes,Vec U,Mat J,Mat B,void *ctx)  */
>>> ierr =
>>> DMCompositeGetLocalVectors(user->packer,&vphi1,&vphi2,&vlambda);CHKERRQ(ierr);
>>>   ierr =
>>> DMCompositeScatter(user->packer,U,vphi1,vphi2,vlambda);CHKERRQ(ierr);
>>>
>>>   ierr = VecGetArray(vlambda,&lambda);CHKERRQ(ierr);
>>>   ierr = DMDAVecGetArray(user->da1,vphi1,&phi1);CHKERRQ(ierr);
>>>   ierr = DMDAVecGetArray(user->da2,vphi2,&phi2);CHKERRQ(ierr);
>>>
>>>   ierr = DMCompositeGetLocalISs(user->packer,&is);CHKERRQ(ierr);
>>>   ierr = MatGetLocalSubMatrix(B,is[0],is[0],&B11);CHKERRQ(ierr);
>>>   ierr = MatGetLocalSubMatrix(B,is[1],is[1],&B22);CHKERRQ(ierr);
>>>
>>>   for (j=ys; j<ys+ym; j++) {
>>>     for (i=xs; i<xs+xm; i++) {
>>>       row.j = j; row.i = i;
>>>       unit = 1.0;
>>>       ierr =
>>> MatSetValuesStencil(B11,1,&row,1,&row,&unit,INSERT_VALUES);CHKERRQ(ierr);
>>>     }
>>>   }
>>>
>>>   for (j=ys; j<ys+ym; j++) {
>>>     for (i=xs; i<xs+xm; i++) {
>>>       row.j = j; row.i = i;
>>>       unit = 1.0;
>>>       ierr =
>>> MatSetValuesStencil(B22,1,&row,1,&row,&unit,INSERT_VALUES);CHKERRQ(ierr);
>>>     }
>>>   }
>>>
>>>   ierr = MatRestoreLocalSubMatrix(B,is[0],is[0],&B11);CHKERRQ(ierr);
>>>   ierr = MatRestoreLocalSubMatrix(B,is[1],is[1],&B22);CHKERRQ(ierr);
>>>
>>>
>>>   unit = 1.0;
>>>   row1  = 960;//last row global index
>>>
>>>   ierr =
>>> MatSetValues(B,1,&row1,1,&row1,&unit,INSERT_VALUES);CHKERRQ(ierr);
>>>
>>>   ierr = VecRestoreArray(vlambda,&lambda);CHKERRQ(ierr);
>>>   ierr = DMDAVecRestoreArray(user->da1,vphi1,&phi1);CHKERRQ(ierr);
>>>   ierr = DMDAVecRestoreArray(user->da2,vphi2,&phi2);CHKERRQ(ierr);
>>> ......
>>>
>>>
>>> Thanks,
>>> Yingjie
>>>
>>