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. Matt > Both versions converge in 3 iterations for the first SNES iteration, > but the relative one starts to churn after that since the residual > starts very small. The true residual goes down to 4/3 and stagnates. > Is there a convenient way to print out the RHS to see whether it has a > component in the nullspace (which seems likely given the true residual > stagnation)? > > I suppose I already do print the result of SNESComputeFunction on the > zero vec, which is > > RHS = [ 6.66666667e-01 1.33333333e+00 6.66666667e-01 > 1.11022302e-16] > > Thanks, > 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
