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

[(Rebecca) Xuefei YUAN]

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