<div dir="ltr"><div dir="ltr"><br><br><div class="gmail_quote"><div dir="ltr">On Mon, Nov 5, 2018 at 10:37 AM Yingjie Wu <<a href="mailto:yjwu16@gmail.com">yjwu16@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><p style="margin:0px;white-space:pre-wrap">Thank you very much for your reply.</p></div><div dir="ltr">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. </div></div></div></div></blockquote><div><br></div><div>OK, DMComposite might be your best choice.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr">I am not very familiar with FieldSplit. I can understand it first.</div><div dir="ltr">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. </div><div dir="ltr">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). <br></div></div></div></div></blockquote><div><br></div><div>So you are trying to create the identity matrix with MatSetValuesStencil but you are getting zeros on the diagonal. </div><div><br></div><div>I would debug this by looking at the matrix in matlab (or Octave). You can use code like this:</div><div><br></div><div><div> if (PETSC_TRUE) {</div><div> PetscViewer viewer;</div><div> ierr = PetscViewerASCIIOpen(comm, "Amat.m", &viewer);CHKERRQ(ierr);</div><div> ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);</div><div> ierr = MatView(Amat,viewer);CHKERRQ(ierr);</div><div> ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);</div><div> ierr = PetscViewerDestroy(&viewer);</div><div> }</div></div><div><br></div><div>Test on a small, serial problem and see where your data is actually going, if not on the diagonal. </div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"></div><div dir="ltr"><br></div><div>Thanks,</div><div>Yingjie</div></div></div></div><br><div class="gmail_quote"><div dir="ltr">Mark Adams <<a href="mailto:mfadams@lbl.gov" target="_blank">mfadams@lbl.gov</a>> 于2018年11月5日周一 下午11:22写道：<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">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.<div><br></div><div>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.</div><div><br></div><div>Mark </div></div><br><div class="gmail_quote"><div dir="ltr">On Mon, Nov 5, 2018 at 10:00 AM Yingjie Wu via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div>Dear Petsc developer:</div><div>Hi,</div><div>I have recently studied the preconditioner of the program, and some problems have arisen. Please help me to solve them. </div><div>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. </div><div>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. </div><div>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. </div></div></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div dir="ltr"><div dir="ltr"><div> yjwu@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</div><div><br></div><div> 0 SNES Function norm 8.235090086536e-02 </div><div>iter = 0, SNES Function norm 0.0823509</div><div>iter = 0, Keff ======= 1.</div><div> Linear solve did not converge due to DIVERGED_PCSETUP_FAILED iterations 0</div><div> PCSETUP_FAILED due to FACTOR_NUMERIC_ZEROPIVOT </div><div>Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0</div><div>SNES Object: 1 MPI processes</div><div> type: newtonls</div><div> maximum iterations=50, maximum function evaluations=10000</div><div> tolerances: relative=1e-08, absolute=1e-50, solution=1e-08</div><div> total number of linear solver iterations=0</div><div> total number of function evaluations=1</div><div> norm schedule ALWAYS</div><div> SNESLineSearch Object: 1 MPI processes</div><div> type: bt</div><div> interpolation: cubic</div><div> alpha=1.000000e-04</div><div> maxstep=1.000000e+08, minlambda=1.000000e-12</div><div> tolerances: relative=1.000000e-08, absolute=1.000000e-15, lambda=1.000000e-08</div><div> maximum iterations=40</div><div> KSP Object: 1 MPI processes</div><div> type: gmres</div><div> restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement</div><div> happy breakdown tolerance 1e-30</div><div> maximum iterations=10000, initial guess is zero</div><div> tolerances: relative=1e-05, absolute=1e-50, divergence=10000.</div><div> left preconditioning</div><div> using PRECONDITIONED norm type for convergence test</div><div> PC Object: 1 MPI processes</div><div> type: ilu</div><div> out-of-place factorization</div><div> 0 levels of fill</div><div> tolerance for zero pivot 2.22045e-14</div><div> matrix ordering: natural</div><div> factor fill ratio given 1., needed 1.</div><div> Factored matrix follows:</div><div> Mat Object: 1 MPI processes</div><div> type: seqaij</div><div> rows=961, cols=961</div><div> package used to perform factorization: petsc</div><div> total: nonzeros=4625, allocated nonzeros=4625</div><div> total number of mallocs used during MatSetValues calls =0</div><div> not using I-node routines</div><div> linear system matrix followed by preconditioner matrix:</div><div> Mat Object: 1 MPI processes</div><div> type: mffd</div><div> rows=961, cols=961</div><div> Matrix-free approximation:</div><div> err=1.49012e-08 (relative error in function evaluation)</div><div> The compute h routine has not yet been set</div><div> Mat Object: 1 MPI processes</div><div> type: seqaij</div><div> rows=961, cols=961</div><div> total: nonzeros=4625, allocated nonzeros=4625</div><div> total number of mallocs used during MatSetValues calls =0</div><div> not using I-node routines</div></div></div></blockquote><div dir="ltr"><div dir="ltr"><div>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. <br></div></div></div></div></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div>......</div></div></div></div></div></blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div>/* set DMComposite */</div></div></div></div></blockquote><div dir="ltr"><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div><div>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);</div><div> ierr = DMSetFromOptions(user.da1);CHKERRQ(ierr);</div><div> ierr = DMSetUp(user.da1);CHKERRQ(ierr);</div><div> ierr = DMCompositeAddDM(user.packer,user.da1);CHKERRQ(ierr);</div><div> 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);</div><div> ierr = DMSetFromOptions(user.da2);CHKERRQ(ierr);</div><div> ierr = DMSetUp(user.da2);CHKERRQ(ierr);</div><div> ierr = DMCompositeAddDM(user.packer,user.da2);CHKERRQ(ierr);</div><div> ierr = DMRedundantCreate(PETSC_COMM_WORLD,0,1,&user.red1);CHKERRQ(ierr);</div><div> ierr = DMCompositeAddDM(user.packer,user.red1);CHKERRQ(ierr);</div></div><div>......</div><div>/* in FormJacobian(SNES snes,Vec U,Mat J,Mat B,void *ctx) */ <br></div><div><div>ierr = DMCompositeGetLocalVectors(user->packer,&vphi1,&vphi2,&vlambda);CHKERRQ(ierr);</div><div> ierr = DMCompositeScatter(user->packer,U,vphi1,vphi2,vlambda);CHKERRQ(ierr);</div><div><br></div><div> ierr = VecGetArray(vlambda,&lambda);CHKERRQ(ierr);</div><div> ierr = DMDAVecGetArray(user->da1,vphi1,&phi1);CHKERRQ(ierr);</div><div> ierr = DMDAVecGetArray(user->da2,vphi2,&phi2);CHKERRQ(ierr);</div><div><br></div><div> ierr = DMCompositeGetLocalISs(user->packer,&is);CHKERRQ(ierr);</div><div> ierr = MatGetLocalSubMatrix(B,is[0],is[0],&B11);CHKERRQ(ierr);</div><div> ierr = MatGetLocalSubMatrix(B,is[1],is[1],&B22);CHKERRQ(ierr);</div><div><br></div></div></div></div></div></blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div><div> for (j=ys; j<ys+ym; j++) {</div><div> for (i=xs; i<xs+xm; i++) {</div><div> row.j = j; row.i = i;</div><div> unit = 1.0;</div><div> ierr = MatSetValuesStencil(B11,1,&row,1,&row,&unit,INSERT_VALUES);CHKERRQ(ierr);</div><div> }</div><div> }</div><div><br></div><div> for (j=ys; j<ys+ym; j++) {</div><div> for (i=xs; i<xs+xm; i++) {</div><div> row.j = j; row.i = i;</div><div> unit = 1.0;</div><div> ierr = MatSetValuesStencil(B22,1,&row,1,&row,&unit,INSERT_VALUES);CHKERRQ(ierr);</div><div> }</div><div> }</div><div><br></div><div> ierr = MatRestoreLocalSubMatrix(B,is[0],is[0],&B11);CHKERRQ(ierr);</div><div> ierr = MatRestoreLocalSubMatrix(B,is[1],is[1],&B22);CHKERRQ(ierr);</div><div><br></div><div><br></div><div> unit = 1.0;</div><div> row1 = 960;//last row global index</div><div><br></div><div> ierr = MatSetValues(B,1,&row1,1,&row1,&unit,INSERT_VALUES);CHKERRQ(ierr);</div><div><br></div><div> ierr = VecRestoreArray(vlambda,&lambda);CHKERRQ(ierr);</div><div> ierr = DMDAVecRestoreArray(user->da1,vphi1,&phi1);CHKERRQ(ierr);</div><div> ierr = DMDAVecRestoreArray(user->da2,vphi2,&phi2);CHKERRQ(ierr);</div></div><div>......</div></div></div></div></blockquote><div dir="ltr"><div dir="ltr"><div dir="ltr"><div><br></div><div>Thanks,</div><div>Yingjie</div></div></div></div></div></div></div></div>
</blockquote></div>
</blockquote></div>
</blockquote></div></div></div>