[petsc-users] difference between solns for np=1 and np>1.

[(Rebecca) Xuefei YUAN]

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