<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Fri, Nov 22, 2013 at 6:42 PM, Geoffrey Irving <span dir="ltr"><<a href="mailto:irving@naml.us" target="_blank">irving@naml.us</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="HOEnZb"><div class="h5">On Fri, Nov 22, 2013 at 4:25 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>
> On Fri, Nov 22, 2013 at 6:09 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>> wrote:<br>
>><br>
>> On Fri, Nov 22, 2013 at 3:41 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
>> wrote:<br>
>> > On Fri, Nov 22, 2013 at 5:36 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>> wrote:<br>
>> >><br>
>> >> I have a duplicate of snes ex12 (FEM Poisson) which works with<br>
>> >> Dirichlet boundary conditions, but it's breaking for me with Neumann<br>
>> >> conditions. In particular, with Neumann conditions I get results<br>
>> >> which explode even though I believe I am setting a constant nullspace.<br>
>> >><br>
>> >> For example, if I use two first order elements (the unit square<br>
>> >> divided into two triangles), the resulting solution has<br>
>> >><br>
>> >> L2 error = 1.75514e+08<br>
>> >> u = [-175513825.75680602, -175513825.66302037,<br>
>> >> -175513825.48390722, -175513824.84436429]<br>
>> >><br>
>> >> This looks rather a lot like the null space isn't getting through. I<br>
>> >> am creating the constant nullspace with<br>
>> >><br>
>> >> MatNullSpace null;<br>
>> >> CHECK(MatNullSpaceCreate(comm(),PETSC_TRUE,0,0,&null));<br>
>> >> CHECK(MatSetNullSpace(m,null));<br>
>> >> CHECK(MatNullSpaceDestroy(&null));<br>
>> >><br>
>> >> If I pass "-ksp_view -mat_view", I get the following. The matrix<br>
>> >> entries seem right (they do indeed have the constant nullspace), and<br>
>> >> ksp_view shows that a nullspace is attached. Is attaching the<br>
>> >> nullspace to the matrix with MatSetNullSpace enough, or do I need to<br>
>> >> additionally attach it to the KSP object?<br>
>> ><br>
>> ><br>
>> > 1) I always run with -ksp_monitor_true_residual now when debugging. This<br>
>> > can<br>
>> > give<br>
>> > you an idea whether you have a singular PC, which I suspect here.<br>
>> ><br>
>> > 2) Can you try using -pc_type jacobi? I think ILU might go crazy on a<br>
>> > deficient matrix.<br>
>><br>
>> Here are results with -ksp_monitor_true_residual -pc_type none:<br>
>><br>
>> <a href="http://naml.us/random/laplace-rtol.txt" target="_blank">http://naml.us/random/laplace-rtol.txt</a> # with -ksp_rtol 1e-5<br>
>> <a href="http://naml.us/random/laplace-atol.txt" target="_blank">http://naml.us/random/laplace-atol.txt</a> # with -ksp_atol 1e-5<br>
><br>
><br>
> Okay, if you have an inconsistent RHS I do not think that true_residual will work<br>
> since it uses the unprojected b, but the solve should be fine.<br>
<br>
</div></div>I still don't understand why the atol version is able to drift so far<br>
away from zero mean, even after tens of thousands of iterations. If<br>
KSP sees a null space on the matrix, shouldn't it project that null<br>
space out of the *linear system* residual and also out of solution on<br>
each iteration? Even if it is only projecting out of the solution<br>
delta, how can null space errors be accumulating?</blockquote><div><br></div><div>Both the KSP and Mat show that the null space is set, so everything should work fine,</div><div>and at this point its no longer DMPlex that is in control, its standard PETSc.</div>
<div><br></div><div>We have reached the limit of usefu talking. Something is obviously wrong with the code,</div><div>but since this routinely works in PETSc examples. In situations like these I think we need</div><div>to follow the execution in the debugger to see what is wrong..You can look at Vec values</div>
<div>in the debugger using</div><div><br></div><div> (gdb) p ((Vec_Seq*) b-.data)->array[0]@v->map.n</div><div><br></div><div>and I look at DMPlex things with</div><div><br></div><div> (gdb) p ((DM_Plex*) dm->data)->coneSection</div>
<div><br></div><div>etc.</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="im">
>> Is there a standard way to do this in the context of an SNES, where<br>
>> I'm not computing the residual directly myself? Should I write a<br>
>> wrapper around DMPlexComputeResidualFEM and pass the wrapper to<br>
>> DMSNESSetFunctionLocal, or is there a way to tell SNES about the<br>
>> nullspace directly? Is such a projection happening somewhere in snes<br>
>> ex12?<br>
<br>
> SNES will do it automatically. You can just call MatNullSpaceRemove().<br>
<br>
</div>I.e., SNES will call MatNullSpaceRemove automatically, or I should<br>
manually in a wrapper around DMSNESSetFunctionLocal?<br>
<br>
Thanks,<br>
Geoffrey<br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>