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