On Sat, Nov 23, 2013 at 4:43 AM, Matthew Knepley <[email protected]> wrote: > On Fri, Nov 22, 2013 at 6:42 PM, Geoffrey Irving <[email protected]> wrote: >> >> On Fri, Nov 22, 2013 at 4:25 PM, Matthew Knepley <[email protected]> >> wrote: >> > On Fri, Nov 22, 2013 at 6:09 PM, Geoffrey Irving <[email protected]> wrote: >> >> >> >> On Fri, Nov 22, 2013 at 3:41 PM, Matthew Knepley <[email protected]> >> >> wrote: >> >> > On Fri, Nov 22, 2013 at 5:36 PM, Geoffrey Irving <[email protected]> >> >> > 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.
Thanks, I appreciate the help. It looks like there were at least two different problems: 1. The boundary FE I was creating had the same dimension as the interior FE (instead of codimension 1), due to misreading ex12 even though I had correctly refactored it. I added a dimension consistency check to my code, but I can do this in DMPlexComputeResidualFEM as well to catch future user errors. 2. Even after fixing the dimensions, my boundary functions in PetscFEM are getting x values both inside and completely outside the domain. Almost certainly more user error, but hopefully also something I can add a check for in petsc once I localize it. Geoffrey
