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

[Geoffrey Irving]

On Fri, Nov 22, 2013 at 4:25 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Fri, Nov 22, 2013 at 6:09 PM, Geoffrey Irving <irving at naml.us> wrote:
>>
>> 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
>
>
> Okay, if you have an inconsistent RHS I do not think that true_residual will work
> since it uses the unprojected b, but the solve should be fine.

I still don't understand why the atol version is able to drift so far
away from zero mean, even after tens of thousands of iterations.  If
KSP sees a null space on the matrix, shouldn't it project that null
space out of the *linear system* residual and also out of solution on
each iteration?  Even if it is only projecting out of the solution
delta, how can null space errors be accumulating?

>> Is there a standard way to do this in the context of an SNES, where
>> I'm not computing the residual directly myself?  Should I write a
>> wrapper around DMPlexComputeResidualFEM and pass the wrapper to
>> DMSNESSetFunctionLocal, or is there a way to tell SNES about the
>> nullspace directly?  Is such a projection happening somewhere in snes
>> ex12?

> SNES will do it automatically. You can just call MatNullSpaceRemove().

I.e., SNES will call MatNullSpaceRemove automatically, or I should
manually in a wrapper around DMSNESSetFunctionLocal?

Thanks,
Geoffrey