Matt, The code fragments they sent sure look like they are using DMDA 2d and they talk about finite differences.
Barry I am sure there are bugs in the unstructured grids code also :-) On May 13, 2014, at 8:27 PM, Matthew Knepley <[email protected]> wrote: > On Tue, May 13, 2014 at 7:55 PM, Barry Smith <[email protected]> wrote: > > What do you mean by ‘''the default “coloring” method’’’ > > If you are using DMDA and either DMGetColoring or the SNESSetDM approach > and dof is 4 then we color each of the 4 variables per grid point with a > different color so coupling between variables within a grid point is not a > problem. This would not explain the problem you are seeing below. > > Run your code with -snes_type test and read the results and follow the > directions to debug your Jacobian. > > I think there may actually be a bug with the coloring for unstructured grids. > I am distilling it down to a nice test case. > > Matt > > > Barry > > > On May 13, 2014, at 1:20 PM, Jean-Arthur Louis Olive <[email protected]> wrote: > > > Hi all, > > we are using PETSc to solve the steady state Stokes equations with > > non-linear viscosities using finite difference. Recently we have realized > > that our true residual norm after the last KSP solve did not match next > > SNES function norm when solving the linear Stokes equations. > > > > So to understand this better, we set up two extremely simple linear > > residuals, one with no coupling between variables (vx, vy, P and T), the > > other with one coupling term (shown below). > > > > RESIDUAL 1 (NO COUPLING): > > for (j=info->ys; j<info->ys+info->ym; j++) { > > for (i=info->xs; i<info->xs+info->xm; i++) { > > f[j][i].P = x[j][i].P - 3000000; > > f[j][i].vx= 2*x[j][i].vx; > > f[j][i].vy= 3*x[j][i].vy - 2; > > f[j][i].T = x[j][i].T; > > } > > > > RESIDUAL 2 (ONE COUPLING TERM): > > for (j=info->ys; j<info->ys+info->ym; j++) { > > for (i=info->xs; i<info->xs+info->xm; i++) { > > f[j][i].P = x[j][i].P - 3; > > f[j][i].vx= x[j][i].vx - 3*x[j][i].vy; > > f[j][i].vy= x[j][i].vy - 2; > > f[j][i].T = x[j][i].T; > > } > > } > > > > > > and our default set of options is: > > > > > > OPTIONS: mpiexec -np $np ../Stokes -snes_max_it 4 -ksp_atol 2.0e+2 > > -ksp_max_it 20 -ksp_rtol 9.0e-1 -ksp_type fgmres -snes_monitor > > -snes_converged_reason -snes_view -log_summary -options_left 1 > > -ksp_monitor_true_residual -pc_type none -snes_linesearch_type cp > > > > > > With the uncoupled residual (Residual 1), we get matching KSP and SNES > > norm, highlighted below: > > > > > > Result from Solve - RESIDUAL 1 > > 0 SNES Function norm 8.485281374240e+07 > > 0 KSP unpreconditioned resid norm 8.485281374240e+07 true resid norm > > 8.485281374240e+07 ||r(i)||/||b|| 1.000000000000e+00 > > 1 KSP unpreconditioned resid norm 1.131370849896e+02 true resid norm > > 1.131370849896e+02 ||r(i)||/||b|| 1.333333333330e-06 > > 1 SNES Function norm 1.131370849896e+02 > > 0 KSP unpreconditioned resid norm 1.131370849896e+02 true resid norm > > 1.131370849896e+02 ||r(i)||/||b|| 1.000000000000e+00 > > 2 SNES Function norm 1.131370849896e+02 > > Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 2 > > > > > > With the coupled residual (Residual 2), the norms do not match, see below > > > > > > Result from Solve - RESIDUAL 2: > > 0 SNES Function norm 1.019803902719e+02 > > 0 KSP unpreconditioned resid norm 1.019803902719e+02 true resid norm > > 1.019803902719e+02 ||r(i)||/||b|| 1.000000000000e+00 > > 1 KSP unpreconditioned resid norm 8.741176309016e+01 true resid norm > > 8.741176309016e+01 ||r(i)||/||b|| 8.571428571429e-01 > > 1 SNES Function norm 1.697056274848e+02 > > 0 KSP unpreconditioned resid norm 1.697056274848e+02 true resid norm > > 1.697056274848e+02 ||r(i)||/||b|| 1.000000000000e+00 > > 1 KSP unpreconditioned resid norm 5.828670868165e-12 true resid norm > > 5.777940247956e-12 ||r(i)||/||b|| 3.404683942184e-14 > > 2 SNES Function norm 3.236770473841e-07 > > Nonlinear solve converged due to CONVERGED_FNORM_RELATIVE iterations 2 > > > > > > Lastly, if we add -snes_fd to our options, the norms for residual 2 get > > better - they match after the first iteration but not after the second. > > > > > > Result from Solve with -snes_fd - RESIDUAL 2 > > 0 SNES Function norm 8.485281374240e+07 > > 0 KSP unpreconditioned resid norm 8.485281374240e+07 true resid norm > > 8.485281374240e+07 ||r(i)||/||b|| 1.000000000000e+00 > > 1 KSP unpreconditioned resid norm 2.039607805429e+02 true resid norm > > 2.039607805429e+02 ||r(i)||/||b|| 2.403700850300e-06 > > 1 SNES Function norm 2.039607805429e+02 > > 0 KSP unpreconditioned resid norm 2.039607805429e+02 true resid norm > > 2.039607805429e+02 ||r(i)||/||b|| 1.000000000000e+00 > > 1 KSP unpreconditioned resid norm 2.529822128436e+01 true resid norm > > 2.529822128436e+01 ||r(i)||/||b|| 1.240347346045e-01 > > 2 SNES Function norm 2.549509757105e+01 [SLIGHTLY DIFFERENT] > > 0 KSP unpreconditioned resid norm 2.549509757105e+01 true resid norm > > 2.549509757105e+01 ||r(i)||/||b|| 1.000000000000e+00 > > 3 SNES Function norm 2.549509757105e+01 > > Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 3 > > > > > > Does this mean that our Jacobian is not approximated properly by the > > default “coloring” method when it has off-diagonal terms? > > > > Thanks a lot, > > Arthur and Eric > > > > > -- > 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
