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) 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 /* 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 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 /* 295 SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created 296 with SNESCreate_NEWTONLS(). 297 298 Input Parameter: 299 . snes - the SNES context 300 301 Application Interface Routine: SNESDestroy() 302 */ 303 static PetscErrorCode SNESDestroy_NEWTONLS(SNES snes) 304 { 305 PetscFunctionBegin; 306 PetscCall(PetscFree(snes->data)); 307 PetscFunctionReturn(PETSC_SUCCESS); 308 } 309 310 /*MC 311 SNESNEWTONLS - Newton based nonlinear solver that uses a line search 312 313 Options Database Keys: 314 + -snes_linesearch_type <bt> - basic (or equivalently none), bt, l2, cp, nleqerr, shell. Select line search type, see `SNESLineSearchSetType()` 315 . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt, see `SNESLineSearchSetOrder()` 316 . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (`SNESLineSearchSetComputeNorms()`) 317 . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient 318 . -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) 319 . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate 320 . -snes_linesearch_monitor - print information about the progress of line searches 321 - -snes_linesearch_damping - damping factor used for basic line search 322 323 Level: beginner 324 325 Note: 326 This is the default nonlinear solver in `SNES` 327 328 .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONTR`, `SNESQN`, `SNESLineSearchSetType()`, `SNESLineSearchSetOrder()` 329 `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` `SNESLineSearchSetComputeNorms()`, `SNESGetLineSearch()`, `SNESLineSearchSetType()` 330 M*/ 331 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes) 332 { 333 SNES_NEWTONLS *neP; 334 SNESLineSearch linesearch; 335 336 PetscFunctionBegin; 337 snes->ops->setup = SNESSetUp_NEWTONLS; 338 snes->ops->solve = SNESSolve_NEWTONLS; 339 snes->ops->destroy = SNESDestroy_NEWTONLS; 340 341 snes->npcside = PC_RIGHT; 342 snes->usesksp = PETSC_TRUE; 343 snes->usesnpc = PETSC_TRUE; 344 345 PetscCall(SNESGetLineSearch(snes, &linesearch)); 346 if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 347 348 snes->alwayscomputesfinalresidual = PETSC_TRUE; 349 350 PetscCall(SNESParametersInitialize(snes)); 351 352 PetscCall(PetscNew(&neP)); 353 snes->data = (void *)neP; 354 PetscFunctionReturn(PETSC_SUCCESS); 355 } 356