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) 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(0); 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, gnorm, xnorm, ynorm; 128 Vec Y, X, F; 129 SNESLineSearch linesearch; 130 SNESConvergedReason reason; 131 132 PetscFunctionBegin; 133 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); 134 135 snes->numFailures = 0; 136 snes->numLinearSolveFailures = 0; 137 snes->reason = SNES_CONVERGED_ITERATING; 138 139 maxits = snes->max_its; /* maximum number of iterations */ 140 X = snes->vec_sol; /* solution vector */ 141 F = snes->vec_func; /* residual vector */ 142 Y = snes->vec_sol_update; /* newton step */ 143 144 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 145 snes->iter = 0; 146 snes->norm = 0.0; 147 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 148 PetscCall(SNESGetLineSearch(snes, &linesearch)); 149 150 /* compute the preconditioned function first in the case of left preconditioning with preconditioned function */ 151 if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) { 152 PetscCall(SNESApplyNPC(snes, X, NULL, F)); 153 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 154 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 155 snes->reason = SNES_DIVERGED_INNER; 156 PetscFunctionReturn(0); 157 } 158 159 PetscCall(VecNormBegin(F, NORM_2, &fnorm)); 160 PetscCall(VecNormEnd(F, NORM_2, &fnorm)); 161 } else { 162 if (!snes->vec_func_init_set) { 163 PetscCall(SNESComputeFunction(snes, X, F)); 164 } else snes->vec_func_init_set = PETSC_FALSE; 165 } 166 167 PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */ 168 SNESCheckFunctionNorm(snes, fnorm); 169 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 170 snes->norm = fnorm; 171 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 172 PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0)); 173 PetscCall(SNESMonitor(snes, 0, fnorm)); 174 175 /* test convergence */ 176 PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP); 177 if (snes->reason) PetscFunctionReturn(0); 178 179 for (i = 0; i < maxits; i++) { 180 /* Call general purpose update function */ 181 PetscTryTypeMethod(snes, update, snes->iter); 182 183 /* apply the nonlinear preconditioner */ 184 if (snes->npc) { 185 if (snes->npcside == PC_RIGHT) { 186 PetscCall(SNESSetInitialFunction(snes->npc, F)); 187 PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0)); 188 PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X)); 189 PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0)); 190 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 191 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 192 snes->reason = SNES_DIVERGED_INNER; 193 PetscFunctionReturn(0); 194 } 195 PetscCall(SNESGetNPCFunction(snes, F, &fnorm)); 196 } else if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) { 197 PetscCall(SNESApplyNPC(snes, X, F, F)); 198 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 199 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 200 snes->reason = SNES_DIVERGED_INNER; 201 PetscFunctionReturn(0); 202 } 203 } 204 } 205 206 /* Solve J Y = F, where J is Jacobian matrix */ 207 PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre)); 208 SNESCheckJacobianDomainerror(snes); 209 PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre)); 210 PetscCall(KSPSolve(snes->ksp, F, Y)); 211 SNESCheckKSPSolve(snes); 212 PetscCall(KSPGetIterationNumber(snes->ksp, &lits)); 213 PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits)); 214 215 if (PetscLogPrintInfo) PetscCall(SNESNEWTONLSCheckResidual_Private(snes, snes->jacobian, F, Y)); 216 217 /* Compute a (scaled) negative update in the line search routine: 218 X <- X - lambda*Y 219 and evaluate F = function(X) (depends on the line search). 220 */ 221 gnorm = fnorm; 222 PetscCall(SNESLineSearchApply(linesearch, X, F, &fnorm, Y)); 223 PetscCall(SNESLineSearchGetReason(linesearch, &lssucceed)); 224 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm)); 225 PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)gnorm, (double)fnorm, (double)ynorm, (int)lssucceed)); 226 if (snes->reason) break; 227 SNESCheckFunctionNorm(snes, fnorm); 228 if (lssucceed) { 229 if (snes->stol * xnorm > ynorm) { 230 snes->reason = SNES_CONVERGED_SNORM_RELATIVE; 231 PetscFunctionReturn(0); 232 } 233 if (++snes->numFailures >= snes->maxFailures) { 234 PetscBool ismin; 235 snes->reason = SNES_DIVERGED_LINE_SEARCH; 236 PetscCall(SNESNEWTONLSCheckLocalMin_Private(snes, snes->jacobian, F, fnorm, &ismin)); 237 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 238 if (snes->errorifnotconverged && snes->reason) { 239 PetscViewer monitor; 240 PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 241 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]); 242 SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due %s.", SNESConvergedReasons[snes->reason]); 243 } 244 break; 245 } 246 } 247 /* Monitor convergence */ 248 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 249 snes->iter = i + 1; 250 snes->norm = fnorm; 251 snes->ynorm = ynorm; 252 snes->xnorm = xnorm; 253 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 254 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits)); 255 PetscCall(SNESMonitor(snes, snes->iter, snes->norm)); 256 /* Test for convergence */ 257 PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP); 258 if (snes->reason) break; 259 } 260 if (i == maxits) { 261 PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits)); 262 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 263 } 264 PetscFunctionReturn(0); 265 } 266 267 /* 268 SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use 269 of the SNESNEWTONLS nonlinear solver. 270 271 Input Parameter: 272 . snes - the SNES context 273 . x - the solution vector 274 275 Application Interface Routine: SNESSetUp() 276 277 */ 278 PetscErrorCode SNESSetUp_NEWTONLS(SNES snes) 279 { 280 PetscFunctionBegin; 281 PetscCall(SNESSetUpMatrices(snes)); 282 if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED; 283 PetscFunctionReturn(0); 284 } 285 286 PetscErrorCode SNESReset_NEWTONLS(SNES snes) 287 { 288 PetscFunctionBegin; 289 PetscFunctionReturn(0); 290 } 291 292 /* 293 SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created 294 with SNESCreate_NEWTONLS(). 295 296 Input Parameter: 297 . snes - the SNES context 298 299 Application Interface Routine: SNESDestroy() 300 */ 301 PetscErrorCode SNESDestroy_NEWTONLS(SNES snes) 302 { 303 PetscFunctionBegin; 304 PetscCall(SNESReset_NEWTONLS(snes)); 305 PetscCall(PetscFree(snes->data)); 306 PetscFunctionReturn(0); 307 } 308 309 /* 310 SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure. 311 312 Input Parameters: 313 . SNES - the SNES context 314 . viewer - visualization context 315 316 Application Interface Routine: SNESView() 317 */ 318 static PetscErrorCode SNESView_NEWTONLS(SNES snes, PetscViewer viewer) 319 { 320 PetscBool iascii; 321 322 PetscFunctionBegin; 323 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 324 if (iascii) { } 325 PetscFunctionReturn(0); 326 } 327 328 /* 329 SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method. 330 331 Input Parameter: 332 . snes - the SNES context 333 334 Application Interface Routine: SNESSetFromOptions() 335 */ 336 static PetscErrorCode SNESSetFromOptions_NEWTONLS(SNES snes, PetscOptionItems *PetscOptionsObject) 337 { 338 PetscFunctionBegin; 339 PetscFunctionReturn(0); 340 } 341 342 /*MC 343 SNESNEWTONLS - Newton based nonlinear solver that uses a line search 344 345 Options Database: 346 + -snes_linesearch_type <bt> - bt,basic. Select line search type 347 . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt 348 . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (`SNESLineSearchSetComputeNorms()`) 349 . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient 350 . -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) 351 . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate 352 . -snes_linesearch_monitor - print information about progress of line searches 353 - -snes_linesearch_damping - damping factor used for basic line search 354 355 Note: 356 This is the default nonlinear solver in `SNES` 357 358 Level: beginner 359 360 .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONTR`, `SNESQN`, `SNESLineSearchSetType()`, `SNESLineSearchSetOrder()` 361 `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` `SNESLineSearchSetComputeNorms()`, `SNESGetLineSearch()` 362 363 M*/ 364 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes) 365 { 366 SNES_NEWTONLS *neP; 367 SNESLineSearch linesearch; 368 369 PetscFunctionBegin; 370 snes->ops->setup = SNESSetUp_NEWTONLS; 371 snes->ops->solve = SNESSolve_NEWTONLS; 372 snes->ops->destroy = SNESDestroy_NEWTONLS; 373 snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS; 374 snes->ops->view = SNESView_NEWTONLS; 375 snes->ops->reset = SNESReset_NEWTONLS; 376 377 snes->npcside = PC_RIGHT; 378 snes->usesksp = PETSC_TRUE; 379 snes->usesnpc = PETSC_TRUE; 380 381 PetscCall(SNESGetLineSearch(snes, &linesearch)); 382 if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 383 384 snes->alwayscomputesfinalresidual = PETSC_TRUE; 385 386 PetscCall(PetscNew(&neP)); 387 snes->data = (void *)neP; 388 PetscFunctionReturn(0); 389 } 390