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

[Matthew Knepley]

On Fri, Nov 22, 2013 at 6:42 PM, Geoffrey Irving <irving at naml.us> wrote:

> 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?


Both the KSP and Mat show that the null space is set, so everything should
work fine,
and at this point its no longer DMPlex that is in control, its standard
PETSc.

We have reached the limit of usefu talking. Something is obviously wrong
with the code,
but since this routinely works in PETSc examples. In situations like these
I think we need
to follow the execution in the debugger to see what is wrong..You can look
at Vec values
in the debugger using

  (gdb) p ((Vec_Seq*) b-.data)->array[0]@v->map.n

and I look at DMPlex things with

  (gdb) p ((DM_Plex*) dm->data)->coneSection

etc.

  Matt


> >> 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
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131123/59d3b1b1/attachment-0001.html>