1 #include <../src/snes/impls/richardson/snesrichardsonimpl.h> 2 3 static PetscErrorCode SNESDestroy_NRichardson(SNES snes) 4 { 5 PetscFunctionBegin; 6 PetscCall(PetscFree(snes->data)); 7 PetscFunctionReturn(PETSC_SUCCESS); 8 } 9 10 static PetscErrorCode SNESSetUp_NRichardson(SNES snes) 11 { 12 PetscFunctionBegin; 13 PetscCheck(snes->npcside != PC_RIGHT, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "NRichardson only supports left preconditioning"); 14 if (snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED; 15 PetscFunctionReturn(PETSC_SUCCESS); 16 } 17 18 static PetscErrorCode SNESSetFromOptions_NRichardson(SNES snes, PetscOptionItems PetscOptionsObject) 19 { 20 PetscFunctionBegin; 21 PetscOptionsHeadBegin(PetscOptionsObject, "SNES Richardson options"); 22 PetscOptionsHeadEnd(); 23 PetscFunctionReturn(PETSC_SUCCESS); 24 } 25 26 static PetscErrorCode SNESSolve_NRichardson(SNES snes) 27 { 28 Vec X, Y, F; 29 PetscReal xnorm, fnorm, ynorm; 30 PetscInt maxits, i; 31 SNESLineSearchReason lsresult; 32 SNESConvergedReason reason; 33 34 PetscFunctionBegin; 35 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); 36 37 snes->reason = SNES_CONVERGED_ITERATING; 38 39 maxits = snes->max_its; /* maximum number of iterations */ 40 X = snes->vec_sol; /* X^n */ 41 Y = snes->vec_sol_update; /* \tilde X */ 42 F = snes->vec_func; /* residual vector */ 43 44 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 45 snes->iter = 0; 46 snes->norm = 0.; 47 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 48 49 if (snes->npc && snes->functype == SNES_FUNCTION_PRECONDITIONED) { 50 PetscCall(SNESApplyNPC(snes, X, NULL, F)); 51 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 52 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 53 snes->reason = SNES_DIVERGED_INNER; 54 PetscFunctionReturn(PETSC_SUCCESS); 55 } 56 PetscCall(VecNorm(F, NORM_2, &fnorm)); 57 } else { 58 if (!snes->vec_func_init_set) PetscCall(SNESComputeFunction(snes, X, F)); 59 else snes->vec_func_init_set = PETSC_FALSE; 60 61 PetscCall(VecNorm(F, NORM_2, &fnorm)); 62 SNESCheckFunctionNorm(snes, fnorm); 63 } 64 if (snes->npc && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) { 65 PetscCall(SNESApplyNPC(snes, X, F, Y)); 66 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 67 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 68 snes->reason = SNES_DIVERGED_INNER; 69 PetscFunctionReturn(PETSC_SUCCESS); 70 } 71 } else { 72 PetscCall(VecCopy(F, Y)); 73 } 74 75 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 76 snes->norm = fnorm; 77 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 78 PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0)); 79 80 /* test convergence */ 81 PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm)); 82 PetscCall(SNESMonitor(snes, 0, fnorm)); 83 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS); 84 85 /* Call general purpose update function */ 86 PetscTryTypeMethod(snes, update, snes->iter); 87 88 for (i = 1; i < maxits + 1; i++) { 89 PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y)); 90 PetscCall(SNESLineSearchGetReason(snes->linesearch, &lsresult)); 91 PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm)); 92 if (lsresult) { 93 if (++snes->numFailures >= snes->maxFailures) { 94 snes->reason = SNES_DIVERGED_LINE_SEARCH; 95 break; 96 } 97 } 98 if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 99 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 100 break; 101 } 102 103 /* Monitor convergence */ 104 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 105 snes->iter = i; 106 snes->norm = fnorm; 107 snes->xnorm = xnorm; 108 snes->ynorm = ynorm; 109 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 110 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0)); 111 /* Test for convergence */ 112 PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm)); 113 PetscCall(SNESMonitor(snes, snes->iter, snes->norm)); 114 if (snes->reason) break; 115 116 /* Call general purpose update function */ 117 PetscTryTypeMethod(snes, update, snes->iter); 118 119 if (snes->npc) { 120 if (snes->functype == SNES_FUNCTION_PRECONDITIONED) { 121 PetscCall(SNESApplyNPC(snes, X, NULL, Y)); 122 PetscCall(VecNorm(F, NORM_2, &fnorm)); 123 PetscCall(VecCopy(Y, F)); 124 } else { 125 PetscCall(SNESApplyNPC(snes, X, F, Y)); 126 } 127 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 128 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 129 snes->reason = SNES_DIVERGED_INNER; 130 PetscFunctionReturn(PETSC_SUCCESS); 131 } 132 } else { 133 PetscCall(VecCopy(F, Y)); 134 } 135 } 136 if (i == maxits + 1) { 137 PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits)); 138 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 139 } 140 PetscFunctionReturn(PETSC_SUCCESS); 141 } 142 143 /*MC 144 SNESNRICHARDSON - Richardson nonlinear solver that uses successive substitutions, also sometimes known as Picard iteration. 145 146 Options Database Keys: 147 + -snes_linesearch_type <l2,cp,basic> - Line search type. 148 - -snes_linesearch_damping <1.0> - Damping for the line search. 149 150 Level: beginner 151 152 Notes: 153 If no inner nonlinear preconditioner is provided then solves $F(x) - b = 0 $ using $x^{n+1} = x^{n} - \lambda 154 (F(x^n) - b) $ where $ \lambda$ is obtained with either `SNESLineSearchSetDamping()`, `-snes_damping` or a line search. If 155 an inner nonlinear preconditioner is provided (either with `-npc_snes_typ`e or `SNESSetNPC()`) then the inner 156 solver is called on the initial solution $x^n$ and the nonlinear Richardson uses $ x^{n+1} = x^{n} + \lambda d^{n}$ 157 where $d^{n} = \hat{x}^{n} - x^{n} $ where $\hat{x}^{n} $ is the solution returned from the inner solver. 158 159 The update, especially without inner nonlinear preconditioner, may be ill-scaled. If using the basic 160 linesearch, one may have to scale the update with `-snes_linesearch_damping` 161 162 This uses no derivative information provided with `SNESSetJacobian()` thus it will be much slower than Newton's method obtained with `-snes_type ls` 163 164 Only supports left non-linear preconditioning. 165 166 .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`, `SNESNGMRES`, `SNESQN`, `SNESNCG`, 167 `SNESLineSearchSetDamping()` 168 M*/ 169 PETSC_EXTERN PetscErrorCode SNESCreate_NRichardson(SNES snes) 170 { 171 SNES_NRichardson *neP; 172 SNESLineSearch linesearch; 173 174 PetscFunctionBegin; 175 snes->ops->destroy = SNESDestroy_NRichardson; 176 snes->ops->setup = SNESSetUp_NRichardson; 177 snes->ops->setfromoptions = SNESSetFromOptions_NRichardson; 178 snes->ops->solve = SNESSolve_NRichardson; 179 180 snes->usesksp = PETSC_FALSE; 181 snes->usesnpc = PETSC_TRUE; 182 183 snes->npcside = PC_LEFT; 184 185 PetscCall(SNESGetLineSearch(snes, &linesearch)); 186 if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHSECANT)); 187 188 snes->alwayscomputesfinalresidual = PETSC_TRUE; 189 190 PetscCall(SNESParametersInitialize(snes)); 191 PetscObjectParameterSetDefault(snes, max_funcs, 30000); 192 PetscObjectParameterSetDefault(snes, max_its, 10000); 193 PetscObjectParameterSetDefault(snes, stol, 1e-20); 194 195 PetscCall(PetscNew(&neP)); 196 snes->data = (void *)neP; 197 PetscFunctionReturn(PETSC_SUCCESS); 198 } 199