[petsc-users] neumann failure in my version of snes ex12

[Geoffrey Irving]

On Fri, Nov 22, 2013 at 3:41 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Fri, Nov 22, 2013 at 5:36 PM, Geoffrey Irving <irving at naml.us> wrote:
>>
>> I have a duplicate of snes ex12 (FEM Poisson) which works with
>> Dirichlet boundary conditions, but it's breaking for me with Neumann
>> conditions.  In particular, with Neumann conditions I get results
>> which explode even though I believe I am setting a constant nullspace.
>>
>> For example, if I use two first order elements (the unit square
>> divided into two triangles), the resulting solution has
>>
>>     L2 error = 1.75514e+08
>>     u = [-175513825.75680602, -175513825.66302037,
>> -175513825.48390722, -175513824.84436429]
>>
>> This looks rather a lot like the null space isn't getting through.  I
>> am creating the constant nullspace with
>>
>>       MatNullSpace null;
>>       CHECK(MatNullSpaceCreate(comm(),PETSC_TRUE,0,0,&null));
>>       CHECK(MatSetNullSpace(m,null));
>>       CHECK(MatNullSpaceDestroy(&null));
>>
>> If I pass "-ksp_view -mat_view", I get the following.  The matrix
>> entries seem right (they do indeed have the constant nullspace), and
>> ksp_view shows that a nullspace is attached.  Is attaching the
>> nullspace to the matrix with MatSetNullSpace enough, or do I need to
>> additionally attach it to the KSP object?
>
>
> 1) I always run with -ksp_monitor_true_residual now when debugging. This can
> give
>     you an idea whether you have a singular PC, which I suspect here.
>
> 2) Can you try using -pc_type jacobi? I think ILU might go crazy on a
> deficient matrix.

Here are results with -ksp_monitor_true_residual -pc_type none:

    http://naml.us/random/laplace-rtol.txt # with -ksp_rtol 1e-5
    http://naml.us/random/laplace-atol.txt # with -ksp_atol 1e-5

Both versions converge in 3 iterations for the first SNES iteration,
but the relative one starts to churn after that since the residual
starts very small.  The true residual goes down to 4/3 and stagnates.
Is there a convenient way to print out the RHS to see whether it has a
component in the nullspace (which seems likely given the true residual
stagnation)?

I suppose I already do print the result of SNESComputeFunction on the
zero vec, which is

    RHS = [  6.66666667e-01   1.33333333e+00   6.66666667e-01   1.11022302e-16]

Thanks,
Geoffrey