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? >> Is there a standard way to do this in the context of an SNES, where >> I'm not computing the residual directly myself? Should I write a >> wrapper around DMPlexComputeResidualFEM and pass the wrapper to >> DMSNESSetFunctionLocal, or is there a way to tell SNES about the >> nullspace directly? Is such a projection happening somewhere in snes >> ex12? > SNES will do it automatically. You can just call MatNullSpaceRemove(). I.e., SNES will call MatNullSpaceRemove automatically, or I should manually in a wrapper around DMSNESSetFunctionLocal? Thanks, Geoffrey
