1 #include <../src/snes/impls/richardson/snesrichardsonimpl.h> 2 3 PetscErrorCode SNESReset_NRichardson(SNES snes) 4 { 5 PetscFunctionBegin; 6 PetscFunctionReturn(PETSC_SUCCESS); 7 } 8 9 /* 10 SNESDestroy_NRichardson - Destroys the private SNES_NRichardson context that was created with SNESCreate_NRichardson(). 11 12 Input Parameter: 13 . snes - the SNES context 14 15 Application Interface Routine: SNESDestroy() 16 */ 17 PetscErrorCode SNESDestroy_NRichardson(SNES snes) 18 { 19 PetscFunctionBegin; 20 PetscCall(SNESReset_NRichardson(snes)); 21 PetscCall(PetscFree(snes->data)); 22 PetscFunctionReturn(PETSC_SUCCESS); 23 } 24 25 /* 26 SNESSetUp_NRichardson - Sets up the internal data structures for the later use 27 of the SNESNRICHARDSON nonlinear solver. 28 29 Input Parameters: 30 + snes - the SNES context 31 - x - the solution vector 32 33 Application Interface Routine: SNESSetUp() 34 */ 35 PetscErrorCode SNESSetUp_NRichardson(SNES snes) 36 { 37 PetscFunctionBegin; 38 PetscCheck(snes->npcside != PC_RIGHT, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "NRichardson only supports left preconditioning"); 39 if (snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED; 40 PetscFunctionReturn(PETSC_SUCCESS); 41 } 42 43 /* 44 SNESSetFromOptions_NRichardson - Sets various parameters for the SNESNEWTONLS method. 45 46 Input Parameter: 47 . snes - the SNES context 48 49 Application Interface Routine: SNESSetFromOptions() 50 */ 51 static PetscErrorCode SNESSetFromOptions_NRichardson(SNES snes, PetscOptionItems *PetscOptionsObject) 52 { 53 PetscFunctionBegin; 54 PetscOptionsHeadBegin(PetscOptionsObject, "SNES Richardson options"); 55 PetscOptionsHeadEnd(); 56 PetscFunctionReturn(PETSC_SUCCESS); 57 } 58 59 /* 60 SNESView_NRichardson - Prints info from the SNESRichardson data structure. 61 62 Input Parameters: 63 + SNES - the SNES context 64 - viewer - visualization context 65 66 Application Interface Routine: SNESView() 67 */ 68 static PetscErrorCode SNESView_NRichardson(SNES snes, PetscViewer viewer) 69 { 70 PetscBool iascii; 71 72 PetscFunctionBegin; 73 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 74 if (iascii) { } 75 PetscFunctionReturn(PETSC_SUCCESS); 76 } 77 78 /* 79 SNESSolve_NRichardson - Solves a nonlinear system with the Richardson method. 80 81 Input Parameters: 82 . snes - the SNES context 83 84 Output Parameter: 85 . outits - number of iterations until termination 86 87 Application Interface Routine: SNESSolve() 88 */ 89 PetscErrorCode SNESSolve_NRichardson(SNES snes) 90 { 91 Vec X, Y, F; 92 PetscReal xnorm, fnorm, ynorm; 93 PetscInt maxits, i; 94 SNESLineSearchReason lsresult; 95 SNESConvergedReason reason; 96 97 PetscFunctionBegin; 98 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); 99 100 snes->reason = SNES_CONVERGED_ITERATING; 101 102 maxits = snes->max_its; /* maximum number of iterations */ 103 X = snes->vec_sol; /* X^n */ 104 Y = snes->vec_sol_update; /* \tilde X */ 105 F = snes->vec_func; /* residual vector */ 106 107 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 108 snes->iter = 0; 109 snes->norm = 0.; 110 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 111 112 if (snes->npc && snes->functype == SNES_FUNCTION_PRECONDITIONED) { 113 PetscCall(SNESApplyNPC(snes, X, NULL, F)); 114 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 115 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 116 snes->reason = SNES_DIVERGED_INNER; 117 PetscFunctionReturn(PETSC_SUCCESS); 118 } 119 PetscCall(VecNorm(F, NORM_2, &fnorm)); 120 } else { 121 if (!snes->vec_func_init_set) { 122 PetscCall(SNESComputeFunction(snes, X, F)); 123 } else snes->vec_func_init_set = PETSC_FALSE; 124 125 PetscCall(VecNorm(F, NORM_2, &fnorm)); 126 SNESCheckFunctionNorm(snes, fnorm); 127 } 128 if (snes->npc && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) { 129 PetscCall(SNESApplyNPC(snes, X, F, Y)); 130 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 131 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 132 snes->reason = SNES_DIVERGED_INNER; 133 PetscFunctionReturn(PETSC_SUCCESS); 134 } 135 } else { 136 PetscCall(VecCopy(F, Y)); 137 } 138 139 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 140 snes->norm = fnorm; 141 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 142 PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0)); 143 PetscCall(SNESMonitor(snes, 0, fnorm)); 144 145 /* test convergence */ 146 PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP); 147 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS); 148 149 /* Call general purpose update function */ 150 PetscTryTypeMethod(snes, update, snes->iter); 151 152 /* set parameter for default relative tolerance convergence test */ 153 snes->ttol = fnorm * snes->rtol; 154 /* test convergence */ 155 PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP); 156 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS); 157 158 for (i = 1; i < maxits + 1; i++) { 159 PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y)); 160 PetscCall(SNESLineSearchGetReason(snes->linesearch, &lsresult)); 161 PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm)); 162 if (lsresult) { 163 if (++snes->numFailures >= snes->maxFailures) { 164 snes->reason = SNES_DIVERGED_LINE_SEARCH; 165 break; 166 } 167 } 168 if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 169 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 170 break; 171 } 172 173 /* Monitor convergence */ 174 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 175 snes->iter = i; 176 snes->norm = fnorm; 177 snes->xnorm = xnorm; 178 snes->ynorm = ynorm; 179 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 180 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0)); 181 PetscCall(SNESMonitor(snes, snes->iter, snes->norm)); 182 /* Test for convergence */ 183 PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP); 184 if (snes->reason) break; 185 186 /* Call general purpose update function */ 187 PetscTryTypeMethod(snes, update, snes->iter); 188 189 if (snes->npc) { 190 if (snes->functype == SNES_FUNCTION_PRECONDITIONED) { 191 PetscCall(SNESApplyNPC(snes, X, NULL, Y)); 192 PetscCall(VecNorm(F, NORM_2, &fnorm)); 193 PetscCall(VecCopy(Y, F)); 194 } else { 195 PetscCall(SNESApplyNPC(snes, X, F, Y)); 196 } 197 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 198 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 199 snes->reason = SNES_DIVERGED_INNER; 200 PetscFunctionReturn(PETSC_SUCCESS); 201 } 202 } else { 203 PetscCall(VecCopy(F, Y)); 204 } 205 } 206 if (i == maxits + 1) { 207 PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits)); 208 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 209 } 210 PetscFunctionReturn(PETSC_SUCCESS); 211 } 212 213 /*MC 214 SNESNRICHARDSON - Richardson nonlinear solver that uses successive substitutions, also sometimes known as Picard iteration. 215 216 Options Database Keys: 217 + -snes_linesearch_type <l2,cp,basic> - Line search type. 218 - -snes_linesearch_damping<1.0> - Damping for the line search. 219 220 Level: beginner 221 222 Notes: 223 If no inner nonlinear preconditioner is provided then solves F(x) - b = 0 using x^{n+1} = x^{n} - lambda 224 (F(x^n) - b) where lambda is obtained either `SNESLineSearchSetDamping()`, -snes_damping or a line search. If 225 an inner nonlinear preconditioner is provided (either with -npc_snes_type or `SNESSetNPC()`) then the inner 226 solver is called an initial solution x^n and the nonlinear Richardson uses x^{n+1} = x^{n} + lambda d^{n} 227 where d^{n} = \hat{x}^{n} - x^{n} where \hat{x}^{n} is the solution returned from the inner solver. 228 229 The update, especially without inner nonlinear preconditioner, may be ill-scaled. If using the basic 230 linesearch, one may have to scale the update with -snes_linesearch_damping 231 232 This uses no derivative information thus will be much slower then Newton's method obtained with -snes_type ls 233 234 Only supports left non-linear preconditioning. 235 236 .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`, `SNESNGMRES`, `SNESQN`, `SNESNCG` 237 M*/ 238 PETSC_EXTERN PetscErrorCode SNESCreate_NRichardson(SNES snes) 239 { 240 SNES_NRichardson *neP; 241 SNESLineSearch linesearch; 242 243 PetscFunctionBegin; 244 snes->ops->destroy = SNESDestroy_NRichardson; 245 snes->ops->setup = SNESSetUp_NRichardson; 246 snes->ops->setfromoptions = SNESSetFromOptions_NRichardson; 247 snes->ops->view = SNESView_NRichardson; 248 snes->ops->solve = SNESSolve_NRichardson; 249 snes->ops->reset = SNESReset_NRichardson; 250 251 snes->usesksp = PETSC_FALSE; 252 snes->usesnpc = PETSC_TRUE; 253 254 snes->npcside = PC_LEFT; 255 256 PetscCall(SNESGetLineSearch(snes, &linesearch)); 257 if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHL2)); 258 259 snes->alwayscomputesfinalresidual = PETSC_TRUE; 260 261 PetscCall(PetscNew(&neP)); 262 snes->data = (void *)neP; 263 264 if (!snes->tolerancesset) { 265 snes->max_funcs = 30000; 266 snes->max_its = 10000; 267 snes->stol = 1e-20; 268 } 269 PetscFunctionReturn(PETSC_SUCCESS); 270 } 271