Dear Jed,
If I checked the Jacobian with
-snes_type test -snes_test_display
np1 has good ratios as:
Norm of matrix ratio 9.88921e-09 difference 3.32649e-06
but
np2 has bad ratios as:
Norm of matrix ratio 1.40616 difference 472.998
I am looking at the hand-coded - finite difference and check those
difference entries.
Let me get back to you later.
BTW, I did look at the tidxm[], tidxn[] and tvv[] in MatSetValues()
ierr = MatSetValues(jac, 1, tidxm, nonzero, tidxn, tvv,
INSERT_VALUES);CHKERRQ(ierr);
in gdb for np1 and np2, they are all the same before calling
ierr =
MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);//book15page37
ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
but why there are quite different in -snes_type test?
Thanks so much for your kind help!
Rebecca
Quoting "(Rebecca) Xuefei YUAN" <xy2102 at columbia.edu>:
> Dear Jed,
>
> I reran the code adding one option -snes_mf, then np=2 gives the same
> result as in np=1.
>
> So I think the bug is in the Jacobian matrix.
>
> Here are steps I used to check on the Jacobian matrix:
>
> 1) After final assemble the Jacobian, output the matrix to Matlab files
> and compare the .m files for both np=1 and np=2:
>
> ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);//book15page37
> ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>
> PetscViewer viewer;
> char fileName[128];
> PetscInt p;
> ierr = MPI_Comm_size(MPI_COMM_WORLD, &p);CHKERRQ(ierr);
> PetscInt its;
> ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
> sprintf(fileName, "twvggt_matrix_tx%i_ty%i_p%i_its%i.m",info1.mx,
> info1.my,p,its);
> ierr =
> PetscViewerASCIIOpen(PETSC_COMM_WORLD,fileName,&viewer);CHKERRQ(ierr);
> ierr =
> PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
> ierr = MatView (jac, viewer); CHKERRQ (ierr);
> PetscViewerDestroy(viewer);
>
> 2) in such a way, after each nonlinear iteration, there is a file for
> the corresponding Jacobian, I found these J_np1 and J_np2 have
> different nonzero structures, for example, for a (7X6+1,7X6+1) matrix,
> standard five-pt scheme with stencil width = 1, (the other 1 after 42
> is a scalar in DMComposite)
> J_np1 has
> % Size = 43 43
> % Nonzeros = 367
> zzz = zeros(367,3);
>
> and J_np2 has
> % Size = 43 43
> % Nonzeros = 576
> zzz = zeros(576,3);
>
> I checked the different nonzero structure in J_np1 and J_np2, the
> different spots in J_np2 are all zeros.
>
> 3) The null space in this case is a single vector has first 7X6
> elements being a scale, and the last element being zero. I confirmed by
> comparing v = null(Mat_0) in Matlab with VecView in C, they are the
> same.
>
> Will this be enough to confirm the J_np1, J_np2 and the null vector?
>
> Thanks very much!
>
> Rebecca
>
>
>
> Quoting Jed Brown <jed at 59A2.org>:
>
>> On Thu, 18 Mar 2010 13:50:32 -0400, "(Rebecca) Xuefei YUAN"
>> <xy2102 at columbia.edu> wrote:
>>> Dear Jed,
>>>
>>> I excluded the bug from CreateNullSpace(), but still have different
>>> convergence history for np=1 and np=2, both with -pc_type none, and
>>> -pc_type jacobi.
>>>
>>> The convergence history for np=1, -pc_type none is
>>>
>>>
>>> 0 SNES Function norm 3.277654936380e+02
>>> Linear solve converged due to CONVERGED_RTOL iterations 1
>>> 1 SNES Function norm 1.010694930474e+01
>>> Linear solve converged due to CONVERGED_RTOL iterations 9
>>> 2 SNES Function norm 1.456202001578e+00
>>> Linear solve converged due to CONVERGED_RTOL iterations 23
>>> 3 SNES Function norm 6.670544108392e-02
>>> Linear solve converged due to CONVERGED_RTOL iterations 28
>>> 4 SNES Function norm 1.924506428876e-04
>>> Linear solve did not converge due to DIVERGED_ITS iterations 10000
>>> 5 SNES Function norm 3.554534723246e-05
>>> Linear solve did not converge due to DIVERGED_ITS iterations 10000
>>> 6 SNES Function norm 3.554534511905e-05
>>> Linear solve did not converge due to DIVERGED_ITS iterations 10000
>>> 7 SNES Function norm 3.554534511895e-05
>>> Nonlinear solve converged due to CONVERGED_PNORM_RELATIVE
>>
>> It looks like this has stagnated. You said you have checked that the
>> matrices are the same, what did you do to confirm this? How did you
>> check that the null spaces are the same? What do the unpreconditioned
>> residuals look like (e.g. -ksp_type fgmres or -ksp_type lgmres
>> -ksp_right_pc)? If you are working in unpreconditioned residuals, then
>> how are you implementing boundary conditions?
>>
>> Jed
>>
>>
>
>
>
> --
> (Rebecca) Xuefei YUAN
> Department of Applied Physics and Applied Mathematics
> Columbia University
> Tel:917-399-8032
> www.columbia.edu/~xy2102
--
(Rebecca) Xuefei YUAN
Department of Applied Physics and Applied Mathematics
Columbia University
Tel:917-399-8032
www.columbia.edu/~xy2102
