<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Sat, Nov 23, 2013 at 12:07 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 Sat, Nov 23, 2013 at 4:43 AM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>
> On Fri, Nov 22, 2013 at 6:42 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>> wrote:<br>
>><br>
>> On Fri, Nov 22, 2013 at 4:25 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
>> 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>><br>
>> >> > 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<br>
>> >> >> 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.<br>
>> >> >> 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.<br>
>> >> > 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<br>
>> > will work<br>
>> > since it uses the unprojected b, but the solve should be fine.<br>
>><br>
>> 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?<br>
><br>
><br>
> Both the KSP and Mat show that the null space is set, so everything should<br>
> work fine,<br>
> and at this point its no longer DMPlex that is in control, its standard<br>
> PETSc.<br>
><br>
> We have reached the limit of usefu talking. Something is obviously wrong with the code,<br>
> but since this routinely works in PETSc examples. In situations like these I think we need<br>
> to follow the execution in the debugger to see what is wrong..You can look at Vec values<br>
> in the debugger using<br>
><br>
> (gdb) p ((Vec_Seq*) b-.data)->array[0]@v->map.n<br>
><br>
> and I look at DMPlex things with<br>
><br>
> (gdb) p ((DM_Plex*) dm->data)->coneSection<br>
><br>
> etc.<br>
<br>
</div></div>Thanks, I appreciate the help. It looks like there were at least two<br>
different problems:<br>
<br>
1. The boundary FE I was creating had the same dimension as the<br>
interior FE (instead of codimension 1), due to misreading ex12 even<br>
though I had correctly refactored it. I added a dimension consistency<br>
check to my code, but I can do this in DMPlexComputeResidualFEM as<br>
well to catch future user errors.<br>
<br>
2. Even after fixing the dimensions, my boundary functions in PetscFEM<br>
are getting x values both inside and completely outside the domain.<br>
Almost certainly more user error, but hopefully also something I can<br>
add a check for in petsc once I localize it.</blockquote><div><br></div><div>This could be my bug. The test I have for ex12 is the variable coefficient problem</div><div>with div (x + y) grad u = f. This seems to check between the analytic and field</div>
<div>versions, meaning that the x coming into f1() matches the x I used to make the</div><div>field, and my exact solution.</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">
<span class="HOEnZb"><font color="#888888"><br>
Geoffrey<br>
</font></span></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>