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
