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