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