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