<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Sun, Nov 24, 2013 at 6:18 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 Sun, Nov 24, 2013 at 4:10 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>
> On Sun, Nov 24, 2013 at 6:06 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>> wrote:<br>
>><br>
>> On Sun, Nov 24, 2013 at 5:41 AM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
>> wrote:<br>
>> > On Sat, Nov 23, 2013 at 5:44 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>> wrote:<br>
>> >><br>
>> >> On Sat, Nov 23, 2013 at 12:20 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
>> >> wrote:<br>
>> >> > On Sat, Nov 23, 2013 at 2:04 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>><br>
>> >> > wrote:<br>
>> >> >><br>
>> >> >> On Sat, Nov 23, 2013 at 10:11 AM, Matthew Knepley<br>
>> >> >> <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
>> >> >> wrote:<br>
>> >> >> > On Sat, Nov 23, 2013 at 12:07 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>><br>
>> >> >> > wrote:<br>
>> >> >> >><br>
>> >> >> >> On Sat, Nov 23, 2013 at 4:43 AM, Matthew Knepley<br>
>> >> >> >> <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
>> >> >> >> wrote:<br>
>> >> >> >> > On Fri, Nov 22, 2013 at 6:42 PM, Geoffrey Irving<br>
>> >> >> >> > <<a href="mailto:irving@naml.us">irving@naml.us</a>><br>
>> >> >> >> > wrote:<br>
>> >> >> >> >><br>
>> >> >> >> >> On Fri, Nov 22, 2013 at 4:25 PM, Matthew Knepley<br>
>> >> >> >> >> <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
>> >> >> >> >> wrote:<br>
>> >> >> >> >> > On Fri, Nov 22, 2013 at 6:09 PM, Geoffrey Irving<br>
>> >> >> >> >> > <<a href="mailto:irving@naml.us">irving@naml.us</a>><br>
>> >> >> >> >> > wrote:<br>
>> >> >> >> >> >><br>
>> >> >> >> >> >> On Fri, Nov 22, 2013 at 3:41 PM, Matthew Knepley<br>
>> >> >> >> >> >> <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
>> >> >> >> >> >> wrote:<br>
>> >> >> >> >> >> > On Fri, Nov 22, 2013 at 5:36 PM, Geoffrey Irving<br>
>> >> >> >> >> >> > <<a href="mailto:irving@naml.us">irving@naml.us</a>><br>
>> >> >> >> >> >> > wrote:<br>
>> >> >> >> >> >> >><br>
>> >> >> >> >> >> >> I have a duplicate of snes ex12 (FEM Poisson) which<br>
>> >> >> >> >> >> >> works<br>
>> >> >> >> >> >> >> with<br>
>> >> >> >> >> >> >> Dirichlet boundary conditions, but it's breaking for me<br>
>> >> >> >> >> >> >> with<br>
>> >> >> >> >> >> >> Neumann<br>
>> >> >> >> >> >> >> conditions. In particular, with Neumann conditions I<br>
>> >> >> >> >> >> >> get<br>
>> >> >> >> >> >> >> results<br>
>> >> >> >> >> >> >> which explode even though I believe I am setting a<br>
>> >> >> >> >> >> >> constant<br>
>> >> >> >> >> >> >> nullspace.<br>
>> >> >> >> >> >> >><br>
>> >> >> >> >> >> >> For example, if I use two first order elements (the unit<br>
>> >> >> >> >> >> >> 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<br>
>> >> >> >> >> >> >> getting<br>
>> >> >> >> >> >> >> through.<br>
>> >> >> >> >> >> >> I<br>
>> >> >> >> >> >> >> am creating the constant nullspace with<br>
>> >> >> >> >> >> >><br>
>> >> >> >> >> >> >> MatNullSpace null;<br>
>> >> >> >> >> >> >><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.<br>
>> >> >> >> >> >> >> The<br>
>> >> >> >> >> >> >> matrix<br>
>> >> >> >> >> >> >> entries seem right (they do indeed have the constant<br>
>> >> >> >> >> >> >> nullspace),<br>
>> >> >> >> >> >> >> and<br>
>> >> >> >> >> >> >> ksp_view shows that a nullspace is attached. Is<br>
>> >> >> >> >> >> >> attaching<br>
>> >> >> >> >> >> >> the<br>
>> >> >> >> >> >> >> nullspace to the matrix with MatSetNullSpace enough, or<br>
>> >> >> >> >> >> >> do<br>
>> >> >> >> >> >> >> I<br>
>> >> >> >> >> >> >> need<br>
>> >> >> >> >> >> >> to<br>
>> >> >> >> >> >> >> additionally attach it to the KSP object?<br>
>> >> >> >> >> >> ><br>
>> >> >> >> >> >> ><br>
>> >> >> >> >> >> > 1) I always run with -ksp_monitor_true_residual now when<br>
>> >> >> >> >> >> > debugging.<br>
>> >> >> >> >> >> > This<br>
>> >> >> >> >> >> > can<br>
>> >> >> >> >> >> > give<br>
>> >> >> >> >> >> > you an idea whether you have a singular PC, which I<br>
>> >> >> >> >> >> > suspect<br>
>> >> >> >> >> >> > here.<br>
>> >> >> >> >> >> ><br>
>> >> >> >> >> >> > 2) Can you try using -pc_type jacobi? I think ILU might<br>
>> >> >> >> >> >> > go<br>
>> >> >> >> >> >> > crazy<br>
>> >> >> >> >> >> > on a<br>
>> >> >> >> >> >> > deficient matrix.<br>
>> >> >> >> >> >><br>
>> >> >> >> >> >> Here are results with -ksp_monitor_true_residual -pc_type<br>
>> >> >> >> >> >> 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<br>
>> >> >> >> >> >> 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<br>
>> >> >> >> >> >> 1e-5<br>
>> >> >> >> >> ><br>
>> >> >> >> >> ><br>
>> >> >> >> >> > Okay, if you have an inconsistent RHS I do not think that<br>
>> >> >> >> >> > true_residual<br>
>> >> >> >> >> > will work<br>
>> >> >> >> >> > since it uses the unprojected b, but the solve should be<br>
>> >> >> >> >> > fine.<br>
>> >> >> >> >><br>
>> >> >> >> >> I still don't understand why the atol version is able to drift<br>
>> >> >> >> >> so<br>
>> >> >> >> >> far<br>
>> >> >> >> >> away from zero mean, even after tens of thousands of<br>
>> >> >> >> >> iterations.<br>
>> >> >> >> >> If<br>
>> >> >> >> >> KSP sees a null space on the matrix, shouldn't it project that<br>
>> >> >> >> >> null<br>
>> >> >> >> >> space out of the *linear system* residual and also out of<br>
>> >> >> >> >> solution<br>
>> >> >> >> >> on<br>
>> >> >> >> >> each iteration? Even if it is only projecting out of the<br>
>> >> >> >> >> 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<br>
>> >> >> >> > everything<br>
>> >> >> >> > should<br>
>> >> >> >> > work fine,<br>
>> >> >> >> > and at this point its no longer DMPlex that is in control, its<br>
>> >> >> >> > standard<br>
>> >> >> >> > PETSc.<br>
>> >> >> >> ><br>
>> >> >> >> > We have reached the limit of usefu talking. Something is<br>
>> >> >> >> > obviously<br>
>> >> >> >> > wrong<br>
>> >> >> >> > with the code,<br>
>> >> >> >> > but since this routinely works in PETSc examples. In situations<br>
>> >> >> >> > like<br>
>> >> >> >> > these I think we need<br>
>> >> >> >> > to follow the execution in the debugger to see what is<br>
>> >> >> >> > wrong..You<br>
>> >> >> >> > can<br>
>> >> >> >> > 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>
>> >> >> >> Thanks, I appreciate the help. It looks like there were at least<br>
>> >> >> >> 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<br>
>> >> >> >> even<br>
>> >> >> >> though I had correctly refactored it. I added a dimension<br>
>> >> >> >> consistency<br>
>> >> >> >> check to my code, but I can do this in DMPlexComputeResidualFEM<br>
>> >> >> >> as<br>
>> >> >> >> well to catch future user errors.<br>
>> >> >> >><br>
>> >> >> >> 2. Even after fixing the dimensions, my boundary functions in<br>
>> >> >> >> PetscFEM<br>
>> >> >> >> are getting x values both inside and completely outside the<br>
>> >> >> >> domain.<br>
>> >> >> >> Almost certainly more user error, but hopefully also something I<br>
>> >> >> >> can<br>
>> >> >> >> add a check for in petsc once I localize it.<br>
>> >> >> ><br>
>> >> >> > This could be my bug. The test I have for ex12 is the variable<br>
>> >> >> > coefficient problem<br>
>> >> >> > with div (x + y) grad u = f. This seems to check between the<br>
>> >> >> > analytic<br>
>> >> >> > and field<br>
>> >> >> > versions, meaning that the x coming into f1() matches the x I used<br>
>> >> >> > to<br>
>> >> >> > make the<br>
>> >> >> > field, and my exact solution.<br>
>> >> >><br>
>> >> >> It does seem to happen with stock snes ex12:<br>
>> >> >><br>
>> >> >> branch: irving/assert-ex12-in-box1<br>
>> >> >><br>
>> >> >> % mpiexec -host localhost -n 1<br>
>> >> >> /home/irving/petsc/debug/lib/ex12-obj/ex12 -run_type test<br>
>> >> >> -refinement_limit 0.0 -bc_type neumann -interpolate 1<br>
>> >> >> -petscspace_order 1 -bd_petscspace_order 1 -show_initial<br>
>> >> >> -dm_plex_print_fem 1 -dm_view ::ascii_info_detail<br>
>> >> >> ...<br>
>> >> >> [0]PETSC ERROR: evaluation at point outside unit box: 0 1.25<br>
>> >> >><br>
>> >> >> I'll trace down why this is happening.<br>
>> >> ><br>
>> >> > My first guess is a triangle with backwards edge. This could cause<br>
>> >> > the<br>
>> >> > geometry routines to barf.<br>
>> >><br>
>> >> I don't think it's edge orientation: it breaks (though at different<br>
>> >> points) regardless of whether I orient all the edges clockwise or<br>
>> >> counterclockwise. Also, I would expect bad edge orientation to result<br>
>> >> in bad normals but not to produce bad quadrature locations (nor bad<br>
>> >> residuals as long as the user routines don't depend on normal).<br>
>> >><br>
>> >> Specifically, I think the problem is a sign error in<br>
>> >> DMPlexComputeProjection2Dto1D_Internal. The following patch seems to<br>
>> >> fix the out of bounds evaluation problem.<br>
>> >> DMPlexComputeProjection2Dto1D is computing a matrix which maps from<br>
>> >> the given segment to the canonical segment, and<br>
>> >> DMPlexComputeLineGeometry_Internal expects a map from the canonical<br>
>> >> segment to the given segment.<br>
>> >><br>
>> >> snes ex12 passes with and without this change, presumably because the<br>
>> >> only Neumann test has constants along each box side, and is therefore<br>
>> >> invariant to this error.<br>
>> ><br>
>> ><br>
>> > When I replaced my kludge with code using the normal explicitly, I get<br>
>> > the<br>
>> > same<br>
>> > error as you do. You fix is correct, and I checked it into<br>
>> > knepley/fix-fem-bd-integrate<br>
>> ><br>
>> > <a href="https://bitbucket.org/petsc/petsc/branch/knepley/fix-fem-bd-integrate" target="_blank">https://bitbucket.org/petsc/petsc/branch/knepley/fix-fem-bd-integrate</a><br>
>> ><br>
>> > along with other fixes that I think get Neumann all the way correct in<br>
>> > ex12.<br>
>> ><br>
>> >><br>
>> >> Unfortunately, my Laplace test is also invariant to this error, so<br>
>> >> this bug is unrelated to the earlier problem.<br>
>> ><br>
>> > It could be one of the other fixes I made. Could you run again?<br>
>><br>
>> Much better. All the points and (inward) normals are right, and now<br>
>> my residuals have zero sum as expected. There's still something<br>
>> wrong, since I don't get the exact solution with second order<br>
>> elements, but that's probably an independent problem. I will trace it<br>
>> down. Thanks for the help and fixes!<br>
>><br>
>> Why did you choose inward pointing normals, by the way? I would have<br>
>> though outward pointing normals are the nearly universal convention.<br>
><br>
> That is a bug fixed in that branch. Did you try next?<br>
<br>
</div></div>Oops, I had misread that commit message as "Use inward pointing<br>
normal" rather than "Used inward pointing normal". I get outward<br>
normals if I reorient the outer boundary edges in my DMPlex to be<br>
clockwise, which I guess is what you intended.<br>
<br>
Closer and closer!</blockquote><div><br></div><div>Hmm, no. I use counter-clockwise ordering too. Now I do not understand what</div><div>is going on. Something is wrong.</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>