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