On Sat, Nov 23, 2013 at 12:07 PM, Geoffrey Irving <[email protected]> wrote:
> 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. This could be my bug. The test I have for ex12 is the variable coefficient problem with div (x + y) grad u = f. This seems to check between the analytic and field versions, meaning that the x coming into f1() matches the x I used to make the field, and my exact solution. Matt > > Geoffrey > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
