Yingjie Wu via petsc-users <[email protected]> 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 <[email protected]> 于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 < >> [email protected]> 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@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 >>> >>
