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 lsreason; 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) 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 SNESCheckFunctionDomainError(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 /* hook state vector to BFGS preconditioner */ 185 PetscCall(KSPGetPC(snes->ksp, &pc)); 186 PetscCall(PCLMVMSetUpdateVec(pc, X)); 187 188 for (i = 0; i < maxits; i++) { 189 /* Call general purpose update function */ 190 PetscTryTypeMethod(snes, update, snes->iter); 191 PetscCall(VecNorm(snes->vec_func, NORM_2, &fnorm)); /* no-op unless update() function changed f() */ 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 if (snes->reason) break; 236 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm)); 237 PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lsreason=%d\n", (double)gnorm, (double)fnorm, (double)ynorm, (int)lsreason)); 238 SNESCheckFunctionDomainError(snes, fnorm); 239 PetscCall(SNESLineSearchGetReason(linesearch, &lsreason)); 240 if (lsreason) { 241 if (snes->stol * xnorm > ynorm) { 242 snes->reason = SNES_CONVERGED_SNORM_RELATIVE; 243 break; 244 } else if (lsreason == SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN) { 245 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 246 break; 247 } else if (lsreason == SNES_LINESEARCH_FAILED_NANORINF) { 248 snes->reason = SNES_DIVERGED_FUNCTION_NANORINF; 249 break; 250 } else if (lsreason == SNES_LINESEARCH_FAILED_OBJECTIVE_DOMAIN) { 251 snes->reason = SNES_DIVERGED_OBJECTIVE_DOMAIN; 252 break; 253 } else if (lsreason == SNES_LINESEARCH_FAILED_JACOBIAN_DOMAIN) { 254 snes->reason = SNES_DIVERGED_JACOBIAN_DOMAIN; 255 break; 256 } else if (++snes->numFailures >= snes->maxFailures) { 257 PetscBool ismin; 258 259 snes->reason = SNES_DIVERGED_LINE_SEARCH; 260 PetscCall(SNESNEWTONLSCheckLocalMin_Private(snes, snes->jacobian, F, fnorm, &ismin)); 261 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 262 break; 263 } 264 } 265 /* Monitor convergence */ 266 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 267 snes->iter = i + 1; 268 snes->norm = fnorm; 269 snes->ynorm = ynorm; 270 snes->xnorm = xnorm; 271 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 272 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits)); 273 /* Test for convergence */ 274 PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm)); 275 PetscCall(SNESMonitor(snes, snes->iter, snes->norm)); 276 if (snes->reason) break; 277 } 278 PetscFunctionReturn(PETSC_SUCCESS); 279 } 280 281 /* 282 SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use 283 of the SNESNEWTONLS nonlinear solver. 284 285 Input Parameter: 286 . snes - the SNES context 287 . x - the solution vector 288 289 Application Interface Routine: SNESSetUp() 290 291 */ 292 static PetscErrorCode SNESSetUp_NEWTONLS(SNES snes) 293 { 294 PetscFunctionBegin; 295 PetscCall(SNESSetUpMatrices(snes)); 296 if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED; 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(PetscFree(snes->data)); 313 PetscFunctionReturn(PETSC_SUCCESS); 314 } 315 316 /*MC 317 SNESNEWTONLS - Newton based nonlinear solver that uses a line search 318 319 Options Database Keys: 320 + -snes_linesearch_type <bt> - basic (or equivalently none), bt, l2, cp, nleqerr, shell. Select line search type, see `SNESLineSearchSetType()` 321 . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt, see `SNESLineSearchSetOrder()` 322 . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (`SNESLineSearchSetComputeNorms()`) 323 . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient 324 . -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) 325 . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate 326 . -snes_linesearch_monitor - print information about the progress of line searches 327 - -snes_linesearch_damping - damping factor used for basic line search 328 329 Level: beginner 330 331 Note: 332 This is the default nonlinear solver in `SNES` 333 334 .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONTR`, `SNESQN`, `SNESLineSearchSetType()`, `SNESLineSearchSetOrder()` 335 `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` `SNESLineSearchSetComputeNorms()`, `SNESGetLineSearch()`, `SNESLineSearchSetType()` 336 M*/ 337 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes) 338 { 339 SNES_NEWTONLS *neP; 340 SNESLineSearch linesearch; 341 342 PetscFunctionBegin; 343 snes->ops->setup = SNESSetUp_NEWTONLS; 344 snes->ops->solve = SNESSolve_NEWTONLS; 345 snes->ops->destroy = SNESDestroy_NEWTONLS; 346 347 snes->npcside = PC_RIGHT; 348 snes->usesksp = PETSC_TRUE; 349 snes->usesnpc = PETSC_TRUE; 350 351 PetscCall(SNESGetLineSearch(snes, &linesearch)); 352 if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 353 354 snes->alwayscomputesfinalresidual = PETSC_TRUE; 355 356 PetscCall(SNESParametersInitialize(snes)); 357 358 PetscCall(PetscNew(&neP)); 359 snes->data = (void *)neP; 360 PetscFunctionReturn(PETSC_SUCCESS); 361 } 362