Hi Matt, our underlying problem is the mismatch between KSP ans SNES norms, even when solving a simple linear system, e.g.,
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;
}
}
which should not have any conditioning issue. So I don’t think in this case
it’s an accuracy problem- but something could be wrong with the FD estimation
of our Jacobian (?)
Arthur
On May 22, 2014, at 11:52 AM, Matthew Knepley <[email protected]> wrote:
> On Thu, May 22, 2014 at 12:47 PM, Jean-Arthur Louis Olive <[email protected]>
> wrote:
> Hi Barry,
> sorry about the late reply-
> We indeed use structured grids (DMDA 2d) - but do not ever provide a Jacobian
> for our non-linear stokes problem (instead just rely on petsc's FD
> approximation). I understand "snes_type test" is meant to compare petsc’s
> Jacobian with a user-provided analytical Jacobian.
> Are you saying we should provide an exact Jacobian for our simple linear test
> and see if there’s a problem with the approximate Jacobian?
>
> The Jacobian computed by PETSc uses a finite-difference approximation, and
> thus is only accurate to maybe 1.0e-7
> depending on the conditioning of your system. Are you trying to compare
> things that are more precise than that? You
> can provide an exact Jacobian to get machine accuracy.
>
> Matt
>
> Thanks,
> Arthur & Eric
>
>
>
> > 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.
> >
> > 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
smime.p7s
Description: S/MIME cryptographic signature
