1 2 #include <../src/snes/impls/ls/lsimpl.h> 3 4 /* 5 Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function, 6 || F(u) ||_2 but not a zero, F(u) = 0. In the case when one cannot compute J^T F we use the fact that 7 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More 8 for this trick. One assumes that the probability that W is in the null space of J is very, very small. 9 */ 10 static PetscErrorCode SNESNEWTONLSCheckLocalMin_Private(SNES snes, Mat A, Vec F, PetscReal fnorm, PetscBool *ismin) { 11 PetscReal a1; 12 PetscBool hastranspose; 13 Vec W; 14 15 PetscFunctionBegin; 16 *ismin = PETSC_FALSE; 17 PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose)); 18 PetscCall(VecDuplicate(F, &W)); 19 if (hastranspose) { 20 /* Compute || J^T F|| */ 21 PetscCall(MatMultTranspose(A, F, W)); 22 PetscCall(VecNorm(W, NORM_2, &a1)); 23 PetscCall(PetscInfo(snes, "|| J^T F|| %14.12e near zero implies found a local minimum\n", (double)(a1 / fnorm))); 24 if (a1 / fnorm < 1.e-4) *ismin = PETSC_TRUE; 25 } else { 26 Vec work; 27 PetscScalar result; 28 PetscReal wnorm; 29 30 PetscCall(VecSetRandom(W, NULL)); 31 PetscCall(VecNorm(W, NORM_2, &wnorm)); 32 PetscCall(VecDuplicate(W, &work)); 33 PetscCall(MatMult(A, W, work)); 34 PetscCall(VecDot(F, work, &result)); 35 PetscCall(VecDestroy(&work)); 36 a1 = PetscAbsScalar(result) / (fnorm * wnorm); 37 PetscCall(PetscInfo(snes, "(F^T J random)/(|| F ||*||random|| %14.12e near zero implies found a local minimum\n", (double)a1)); 38 if (a1 < 1.e-4) *ismin = PETSC_TRUE; 39 } 40 PetscCall(VecDestroy(&W)); 41 PetscFunctionReturn(0); 42 } 43 44 /* 45 Checks if J^T(F - J*X) = 0 46 */ 47 static PetscErrorCode SNESNEWTONLSCheckResidual_Private(SNES snes, Mat A, Vec F, Vec X) { 48 PetscReal a1, a2; 49 PetscBool hastranspose; 50 51 PetscFunctionBegin; 52 PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose)); 53 if (hastranspose) { 54 Vec W1, W2; 55 56 PetscCall(VecDuplicate(F, &W1)); 57 PetscCall(VecDuplicate(F, &W2)); 58 PetscCall(MatMult(A, X, W1)); 59 PetscCall(VecAXPY(W1, -1.0, F)); 60 61 /* Compute || J^T W|| */ 62 PetscCall(MatMultTranspose(A, W1, W2)); 63 PetscCall(VecNorm(W1, NORM_2, &a1)); 64 PetscCall(VecNorm(W2, NORM_2, &a2)); 65 if (a1 != 0.0) PetscCall(PetscInfo(snes, "||J^T(F-Ax)||/||F-AX|| %14.12e near zero implies inconsistent rhs\n", (double)(a2 / a1))); 66 PetscCall(VecDestroy(&W1)); 67 PetscCall(VecDestroy(&W2)); 68 } 69 PetscFunctionReturn(0); 70 } 71 72 /* -------------------------------------------------------------------- 73 74 This file implements a truncated Newton method with a line search, 75 for solving a system of nonlinear equations, using the KSP, Vec, 76 and Mat interfaces for linear solvers, vectors, and matrices, 77 respectively. 78 79 The following basic routines are required for each nonlinear solver: 80 SNESCreate_XXX() - Creates a nonlinear solver context 81 SNESSetFromOptions_XXX() - Sets runtime options 82 SNESSolve_XXX() - Solves the nonlinear system 83 SNESDestroy_XXX() - Destroys the nonlinear solver context 84 The suffix "_XXX" denotes a particular implementation, in this case 85 we use _NEWTONLS (e.g., SNESCreate_NEWTONLS, SNESSolve_NEWTONLS) for solving 86 systems of nonlinear equations with a line search (LS) method. 87 These routines are actually called via the common user interface 88 routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and 89 SNESDestroy(), so the application code interface remains identical 90 for all nonlinear solvers. 91 92 Another key routine is: 93 SNESSetUp_XXX() - Prepares for the use of a nonlinear solver 94 by setting data structures and options. The interface routine SNESSetUp() 95 is not usually called directly by the user, but instead is called by 96 SNESSolve() if necessary. 97 98 Additional basic routines are: 99 SNESView_XXX() - Prints details of runtime options that 100 have actually been used. 101 These are called by application codes via the interface routines 102 SNESView(). 103 104 The various types of solvers (preconditioners, Krylov subspace methods, 105 nonlinear solvers, timesteppers) are all organized similarly, so the 106 above description applies to these categories also. 107 108 -------------------------------------------------------------------- */ 109 /* 110 SNESSolve_NEWTONLS - Solves a nonlinear system with a truncated Newton 111 method with a line search. 112 113 Input Parameters: 114 . snes - the SNES context 115 116 Output Parameter: 117 . outits - number of iterations until termination 118 119 Application Interface Routine: SNESSolve() 120 121 Notes: 122 This implements essentially a truncated Newton method with a 123 line search. By default a cubic backtracking line search 124 is employed, as described in the text "Numerical Methods for 125 Unconstrained Optimization and Nonlinear Equations" by Dennis 126 and Schnabel. 127 */ 128 PetscErrorCode SNESSolve_NEWTONLS(SNES snes) { 129 PetscInt maxits, i, lits; 130 SNESLineSearchReason lssucceed; 131 PetscReal fnorm, gnorm, xnorm, ynorm; 132 Vec Y, X, F; 133 SNESLineSearch linesearch; 134 SNESConvergedReason reason; 135 136 PetscFunctionBegin; 137 PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name); 138 139 snes->numFailures = 0; 140 snes->numLinearSolveFailures = 0; 141 snes->reason = SNES_CONVERGED_ITERATING; 142 143 maxits = snes->max_its; /* maximum number of iterations */ 144 X = snes->vec_sol; /* solution vector */ 145 F = snes->vec_func; /* residual vector */ 146 Y = snes->vec_sol_update; /* newton step */ 147 148 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 149 snes->iter = 0; 150 snes->norm = 0.0; 151 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 152 PetscCall(SNESGetLineSearch(snes, &linesearch)); 153 154 /* compute the preconditioned function first in the case of left preconditioning with preconditioned function */ 155 if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) { 156 PetscCall(SNESApplyNPC(snes, X, NULL, F)); 157 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 158 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 159 snes->reason = SNES_DIVERGED_INNER; 160 PetscFunctionReturn(0); 161 } 162 163 PetscCall(VecNormBegin(F, NORM_2, &fnorm)); 164 PetscCall(VecNormEnd(F, NORM_2, &fnorm)); 165 } else { 166 if (!snes->vec_func_init_set) { 167 PetscCall(SNESComputeFunction(snes, X, F)); 168 } else snes->vec_func_init_set = PETSC_FALSE; 169 } 170 171 PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */ 172 SNESCheckFunctionNorm(snes, fnorm); 173 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 174 snes->norm = fnorm; 175 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 176 PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0)); 177 PetscCall(SNESMonitor(snes, 0, fnorm)); 178 179 /* test convergence */ 180 PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP); 181 if (snes->reason) PetscFunctionReturn(0); 182 183 for (i = 0; i < maxits; i++) { 184 /* Call general purpose update function */ 185 PetscTryTypeMethod(snes, update, snes->iter); 186 187 /* apply the nonlinear preconditioner */ 188 if (snes->npc) { 189 if (snes->npcside == PC_RIGHT) { 190 PetscCall(SNESSetInitialFunction(snes->npc, F)); 191 PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0)); 192 PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X)); 193 PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0)); 194 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 195 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 196 snes->reason = SNES_DIVERGED_INNER; 197 PetscFunctionReturn(0); 198 } 199 PetscCall(SNESGetNPCFunction(snes, F, &fnorm)); 200 } else if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) { 201 PetscCall(SNESApplyNPC(snes, X, F, F)); 202 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 203 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 204 snes->reason = SNES_DIVERGED_INNER; 205 PetscFunctionReturn(0); 206 } 207 } 208 } 209 210 /* Solve J Y = F, where J is Jacobian matrix */ 211 PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre)); 212 SNESCheckJacobianDomainerror(snes); 213 PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre)); 214 PetscCall(KSPSolve(snes->ksp, F, Y)); 215 SNESCheckKSPSolve(snes); 216 PetscCall(KSPGetIterationNumber(snes->ksp, &lits)); 217 PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits)); 218 219 if (PetscLogPrintInfo) PetscCall(SNESNEWTONLSCheckResidual_Private(snes, snes->jacobian, F, Y)); 220 221 /* Compute a (scaled) negative update in the line search routine: 222 X <- X - lambda*Y 223 and evaluate F = function(X) (depends on the line search). 224 */ 225 gnorm = fnorm; 226 PetscCall(SNESLineSearchApply(linesearch, X, F, &fnorm, Y)); 227 PetscCall(SNESLineSearchGetReason(linesearch, &lssucceed)); 228 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm)); 229 PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)gnorm, (double)fnorm, (double)ynorm, (int)lssucceed)); 230 if (snes->reason) break; 231 SNESCheckFunctionNorm(snes, fnorm); 232 if (lssucceed) { 233 if (snes->stol * xnorm > ynorm) { 234 snes->reason = SNES_CONVERGED_SNORM_RELATIVE; 235 PetscFunctionReturn(0); 236 } 237 if (++snes->numFailures >= snes->maxFailures) { 238 PetscBool ismin; 239 snes->reason = SNES_DIVERGED_LINE_SEARCH; 240 PetscCall(SNESNEWTONLSCheckLocalMin_Private(snes, snes->jacobian, F, fnorm, &ismin)); 241 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 242 if (snes->errorifnotconverged && snes->reason) { 243 PetscViewer monitor; 244 PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 245 PetscCheck(monitor, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due to %s. Suggest running with -snes_linesearch_monitor", SNESConvergedReasons[snes->reason]); 246 SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due %s.", SNESConvergedReasons[snes->reason]); 247 } 248 break; 249 } 250 } 251 /* Monitor convergence */ 252 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 253 snes->iter = i + 1; 254 snes->norm = fnorm; 255 snes->ynorm = ynorm; 256 snes->xnorm = xnorm; 257 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 258 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits)); 259 PetscCall(SNESMonitor(snes, snes->iter, snes->norm)); 260 /* Test for convergence */ 261 PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP); 262 if (snes->reason) break; 263 } 264 if (i == maxits) { 265 PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits)); 266 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 267 } 268 PetscFunctionReturn(0); 269 } 270 /* -------------------------------------------------------------------------- */ 271 /* 272 SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use 273 of the SNESNEWTONLS nonlinear solver. 274 275 Input Parameter: 276 . snes - the SNES context 277 . x - the solution vector 278 279 Application Interface Routine: SNESSetUp() 280 281 Notes: 282 For basic use of the SNES solvers, the user need not explicitly call 283 SNESSetUp(), since these actions will automatically occur during 284 the call to SNESSolve(). 285 */ 286 PetscErrorCode SNESSetUp_NEWTONLS(SNES snes) { 287 PetscFunctionBegin; 288 PetscCall(SNESSetUpMatrices(snes)); 289 if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED; 290 PetscFunctionReturn(0); 291 } 292 /* -------------------------------------------------------------------------- */ 293 294 PetscErrorCode SNESReset_NEWTONLS(SNES snes) { 295 PetscFunctionBegin; 296 PetscFunctionReturn(0); 297 } 298 299 /* 300 SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created 301 with SNESCreate_NEWTONLS(). 302 303 Input Parameter: 304 . snes - the SNES context 305 306 Application Interface Routine: SNESDestroy() 307 */ 308 PetscErrorCode SNESDestroy_NEWTONLS(SNES snes) { 309 PetscFunctionBegin; 310 PetscCall(SNESReset_NEWTONLS(snes)); 311 PetscCall(PetscFree(snes->data)); 312 PetscFunctionReturn(0); 313 } 314 /* -------------------------------------------------------------------------- */ 315 316 /* 317 SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure. 318 319 Input Parameters: 320 . SNES - the SNES context 321 . viewer - visualization context 322 323 Application Interface Routine: SNESView() 324 */ 325 static PetscErrorCode SNESView_NEWTONLS(SNES snes, PetscViewer viewer) { 326 PetscBool iascii; 327 328 PetscFunctionBegin; 329 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 330 if (iascii) { } 331 PetscFunctionReturn(0); 332 } 333 334 /* -------------------------------------------------------------------------- */ 335 /* 336 SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method. 337 338 Input Parameter: 339 . snes - the SNES context 340 341 Application Interface Routine: SNESSetFromOptions() 342 */ 343 static PetscErrorCode SNESSetFromOptions_NEWTONLS(SNES snes, PetscOptionItems *PetscOptionsObject) { 344 PetscFunctionBegin; 345 PetscFunctionReturn(0); 346 } 347 348 /* -------------------------------------------------------------------------- */ 349 /*MC 350 SNESNEWTONLS - Newton based nonlinear solver that uses a line search 351 352 Options Database: 353 + -snes_linesearch_type <bt> - bt,basic. Select line search type 354 . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt 355 . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (SNESLineSearchSetComputeNorms()) 356 . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient 357 . -snes_linesearch_maxstep <maxstep> - Sets the maximum stepsize the line search will use (if the 2-norm(y) > maxstep then scale y to be y = (maxstep/2-norm(y)) *y) 358 . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate 359 . -snes_linesearch_monitor - print information about progress of line searches 360 - -snes_linesearch_damping - damping factor used for basic line search 361 362 Notes: 363 This is the default nonlinear solver in SNES 364 365 Level: beginner 366 367 .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONTR`, `SNESQN`, `SNESLineSearchSetType()`, `SNESLineSearchSetOrder()` 368 `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` `SNESLineSearchSetComputeNorms()` 369 370 M*/ 371 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes) { 372 SNES_NEWTONLS *neP; 373 SNESLineSearch linesearch; 374 375 PetscFunctionBegin; 376 snes->ops->setup = SNESSetUp_NEWTONLS; 377 snes->ops->solve = SNESSolve_NEWTONLS; 378 snes->ops->destroy = SNESDestroy_NEWTONLS; 379 snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS; 380 snes->ops->view = SNESView_NEWTONLS; 381 snes->ops->reset = SNESReset_NEWTONLS; 382 383 snes->npcside = PC_RIGHT; 384 snes->usesksp = PETSC_TRUE; 385 snes->usesnpc = PETSC_TRUE; 386 387 PetscCall(SNESGetLineSearch(snes, &linesearch)); 388 if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 389 390 snes->alwayscomputesfinalresidual = PETSC_TRUE; 391 392 PetscCall(PetscNewLog(snes, &neP)); 393 snes->data = (void *)neP; 394 PetscFunctionReturn(0); 395 } 396