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 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 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 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 /* set parameter for default relative tolerance convergence test */ 107 snes->ttol = fnorm * snes->rtol; 108 /* test convergence */ 109 PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm)); 110 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS); 111 112 for (i = 1; i < maxits + 1; i++) { 113 PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y)); 114 PetscCall(SNESLineSearchGetReason(snes->linesearch, &lsresult)); 115 PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm)); 116 if (lsresult) { 117 if (++snes->numFailures >= snes->maxFailures) { 118 snes->reason = SNES_DIVERGED_LINE_SEARCH; 119 break; 120 } 121 } 122 if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 123 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 124 break; 125 } 126 127 /* Monitor convergence */ 128 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 129 snes->iter = i; 130 snes->norm = fnorm; 131 snes->xnorm = xnorm; 132 snes->ynorm = ynorm; 133 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 134 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0)); 135 /* Test for convergence */ 136 PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm)); 137 PetscCall(SNESMonitor(snes, snes->iter, snes->norm)); 138 if (snes->reason) break; 139 140 /* Call general purpose update function */ 141 PetscTryTypeMethod(snes, update, snes->iter); 142 143 if (snes->npc) { 144 if (snes->functype == SNES_FUNCTION_PRECONDITIONED) { 145 PetscCall(SNESApplyNPC(snes, X, NULL, Y)); 146 PetscCall(VecNorm(F, NORM_2, &fnorm)); 147 PetscCall(VecCopy(Y, F)); 148 } else { 149 PetscCall(SNESApplyNPC(snes, X, F, Y)); 150 } 151 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 152 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 153 snes->reason = SNES_DIVERGED_INNER; 154 PetscFunctionReturn(PETSC_SUCCESS); 155 } 156 } else { 157 PetscCall(VecCopy(F, Y)); 158 } 159 } 160 if (i == maxits + 1) { 161 PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits)); 162 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 163 } 164 PetscFunctionReturn(PETSC_SUCCESS); 165 } 166 167 /*MC 168 SNESNRICHARDSON - Richardson nonlinear solver that uses successive substitutions, also sometimes known as Picard iteration. 169 170 Options Database Keys: 171 + -snes_linesearch_type <l2,cp,basic> - Line search type. 172 - -snes_linesearch_damping<1.0> - Damping for the line search. 173 174 Level: beginner 175 176 Notes: 177 If no inner nonlinear preconditioner is provided then solves F(x) - b = 0 using x^{n+1} = x^{n} - lambda 178 (F(x^n) - b) where lambda is obtained either `SNESLineSearchSetDamping()`, -snes_damping or a line search. If 179 an inner nonlinear preconditioner is provided (either with -npc_snes_type or `SNESSetNPC()`) then the inner 180 solver is called an initial solution x^n and the nonlinear Richardson uses x^{n+1} = x^{n} + lambda d^{n} 181 where d^{n} = \hat{x}^{n} - x^{n} where \hat{x}^{n} is the solution returned from the inner solver. 182 183 The update, especially without inner nonlinear preconditioner, may be ill-scaled. If using the basic 184 linesearch, one may have to scale the update with -snes_linesearch_damping 185 186 This uses no derivative information thus will be much slower then Newton's method obtained with -snes_type ls 187 188 Only supports left non-linear preconditioning. 189 190 .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`, `SNESNGMRES`, `SNESQN`, `SNESNCG` 191 M*/ 192 PETSC_EXTERN PetscErrorCode SNESCreate_NRichardson(SNES snes) 193 { 194 SNES_NRichardson *neP; 195 SNESLineSearch linesearch; 196 197 PetscFunctionBegin; 198 snes->ops->destroy = SNESDestroy_NRichardson; 199 snes->ops->setup = SNESSetUp_NRichardson; 200 snes->ops->setfromoptions = SNESSetFromOptions_NRichardson; 201 snes->ops->view = SNESView_NRichardson; 202 snes->ops->solve = SNESSolve_NRichardson; 203 snes->ops->reset = SNESReset_NRichardson; 204 205 snes->usesksp = PETSC_FALSE; 206 snes->usesnpc = PETSC_TRUE; 207 208 snes->npcside = PC_LEFT; 209 210 PetscCall(SNESGetLineSearch(snes, &linesearch)); 211 if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHL2)); 212 213 snes->alwayscomputesfinalresidual = PETSC_TRUE; 214 215 PetscCall(PetscNew(&neP)); 216 snes->data = (void *)neP; 217 218 if (!snes->tolerancesset) { 219 snes->max_funcs = 30000; 220 snes->max_its = 10000; 221 snes->stol = 1e-20; 222 } 223 PetscFunctionReturn(PETSC_SUCCESS); 224 } 225