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