<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Fri, Nov 22, 2013 at 6:09 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">On Fri, Nov 22, 2013 at 3:41 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> 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 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></blockquote><div><br></div><div>Okay, if you have an inconsistent RHS I do not think that true_residual will work</div>
<div>since it uses the unprojected b, but the solve should be fine.</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">
Both versions converge in 3 iterations for the first SNES iteration,<br>
but the relative one starts to churn after that since the residual<br>
starts very small.  The true residual goes down to 4/3 and stagnates.<br>
Is there a convenient way to print out the RHS to see whether it has a<br>
component in the nullspace (which seems likely given the true residual<br>
stagnation)?<br>
<br>
I suppose I already do print the result of SNESComputeFunction on the<br>
zero vec, which is<br>
<br>
    RHS = [  6.66666667e-01   1.33333333e+00   6.66666667e-01   1.11022302e-16]<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>