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