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