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