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 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 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 && reason != SNES_DIVERGED_TR_DELTA) { 202 PetscCall(SNESSetConvergedReason(snes, SNES_DIVERGED_INNER)); 203 PetscFunctionReturn(PETSC_SUCCESS); 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 && reason != SNES_DIVERGED_TR_DELTA) { 210 PetscCall(SNESSetConvergedReason(snes, SNES_DIVERGED_INNER)); 211 PetscFunctionReturn(PETSC_SUCCESS); 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 #if defined(PETSC_USE_INFO) 228 gnorm = fnorm; 229 #endif 230 /* Compute a (scaled) negative update in the line search routine: 231 X <- X - lambda*Y 232 and evaluate F = function(X) (depends on the line search). 233 */ 234 PetscCall(SNESLineSearchApply(linesearch, X, F, &fnorm, Y)); 235 PetscCall(SNESLineSearchGetReason(linesearch, &lssucceed)); 236 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm)); 237 PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)gnorm, (double)fnorm, (double)ynorm, (int)lssucceed)); 238 if (snes->reason) break; 239 SNESCheckFunctionNorm(snes, fnorm); 240 if (lssucceed) { 241 if (snes->stol * xnorm > ynorm) { 242 snes->reason = SNES_CONVERGED_SNORM_RELATIVE; 243 PetscFunctionReturn(PETSC_SUCCESS); 244 } 245 if (++snes->numFailures >= snes->maxFailures) { 246 PetscBool ismin; 247 snes->reason = SNES_DIVERGED_LINE_SEARCH; 248 PetscCall(SNESNEWTONLSCheckLocalMin_Private(snes, snes->jacobian, F, fnorm, &ismin)); 249 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 250 if (snes->errorifnotconverged && snes->reason) { 251 PetscViewer monitor; 252 PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 253 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]); 254 SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due %s.", SNESConvergedReasons[snes->reason]); 255 } 256 break; 257 } 258 } 259 /* Monitor convergence */ 260 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 261 snes->iter = i + 1; 262 snes->norm = fnorm; 263 snes->ynorm = ynorm; 264 snes->xnorm = xnorm; 265 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 266 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits)); 267 /* Test for convergence */ 268 PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm)); 269 PetscCall(SNESMonitor(snes, snes->iter, snes->norm)); 270 if (snes->reason) break; 271 } 272 PetscFunctionReturn(PETSC_SUCCESS); 273 } 274 275 /* 276 SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use 277 of the SNESNEWTONLS nonlinear solver. 278 279 Input Parameter: 280 . snes - the SNES context 281 . x - the solution vector 282 283 Application Interface Routine: SNESSetUp() 284 285 */ 286 static PetscErrorCode SNESSetUp_NEWTONLS(SNES snes) 287 { 288 PetscFunctionBegin; 289 PetscCall(SNESSetUpMatrices(snes)); 290 if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED; 291 PetscFunctionReturn(PETSC_SUCCESS); 292 } 293 294 static PetscErrorCode SNESReset_NEWTONLS(SNES snes) 295 { 296 PetscFunctionBegin; 297 PetscFunctionReturn(PETSC_SUCCESS); 298 } 299 300 /* 301 SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created 302 with SNESCreate_NEWTONLS(). 303 304 Input Parameter: 305 . snes - the SNES context 306 307 Application Interface Routine: SNESDestroy() 308 */ 309 static PetscErrorCode SNESDestroy_NEWTONLS(SNES snes) 310 { 311 PetscFunctionBegin; 312 PetscCall(SNESReset_NEWTONLS(snes)); 313 PetscCall(PetscFree(snes->data)); 314 PetscFunctionReturn(PETSC_SUCCESS); 315 } 316 317 /* 318 SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure. 319 320 Input Parameters: 321 . SNES - the SNES context 322 . viewer - visualization context 323 324 Application Interface Routine: SNESView() 325 */ 326 static PetscErrorCode SNESView_NEWTONLS(SNES snes, PetscViewer viewer) 327 { 328 PetscBool iascii; 329 330 PetscFunctionBegin; 331 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 332 if (iascii) { } 333 PetscFunctionReturn(PETSC_SUCCESS); 334 } 335 336 /* 337 SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method. 338 339 Input Parameter: 340 . snes - the SNES context 341 342 Application Interface Routine: SNESSetFromOptions() 343 */ 344 static PetscErrorCode SNESSetFromOptions_NEWTONLS(SNES snes, PetscOptionItems *PetscOptionsObject) 345 { 346 PetscFunctionBegin; 347 PetscFunctionReturn(PETSC_SUCCESS); 348 } 349 350 /*MC 351 SNESNEWTONLS - Newton based nonlinear solver that uses a line search 352 353 Options Database Keys: 354 + -snes_linesearch_type <bt> - bt,basic. Select line search type 355 . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt 356 . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (`SNESLineSearchSetComputeNorms()`) 357 . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient 358 . -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) 359 . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate 360 . -snes_linesearch_monitor - print information about progress of line searches 361 - -snes_linesearch_damping - damping factor used for basic line search 362 363 Level: beginner 364 365 Note: 366 This is the default nonlinear solver in `SNES` 367 368 .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONTR`, `SNESQN`, `SNESLineSearchSetType()`, `SNESLineSearchSetOrder()` 369 `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` `SNESLineSearchSetComputeNorms()`, `SNESGetLineSearch()` 370 M*/ 371 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes) 372 { 373 SNES_NEWTONLS *neP; 374 SNESLineSearch linesearch; 375 376 PetscFunctionBegin; 377 snes->ops->setup = SNESSetUp_NEWTONLS; 378 snes->ops->solve = SNESSolve_NEWTONLS; 379 snes->ops->destroy = SNESDestroy_NEWTONLS; 380 snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS; 381 snes->ops->view = SNESView_NEWTONLS; 382 snes->ops->reset = SNESReset_NEWTONLS; 383 384 snes->npcside = PC_RIGHT; 385 snes->usesksp = PETSC_TRUE; 386 snes->usesnpc = PETSC_TRUE; 387 388 PetscCall(SNESGetLineSearch(snes, &linesearch)); 389 if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 390 391 snes->alwayscomputesfinalresidual = PETSC_TRUE; 392 393 PetscCall(PetscNew(&neP)); 394 snes->data = (void *)neP; 395 PetscFunctionReturn(PETSC_SUCCESS); 396 } 397