1 #include <../src/snes/impls/richardson/snesrichardsonimpl.h> 2 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "SNESReset_NRichardson" 6 PetscErrorCode SNESReset_NRichardson(SNES snes) 7 { 8 SNES_NRichardson *neP = (SNES_NRichardson *)snes->data; 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 ierr = PetscLineSearchDestroy(&neP->linesearch);CHKERRQ(ierr); 13 PetscFunctionReturn(0); 14 } 15 16 /* 17 SNESDestroy_NRichardson - Destroys the private SNES_NRichardson context that was created with SNESCreate_NRichardson(). 18 19 Input Parameter: 20 . snes - the SNES context 21 22 Application Interface Routine: SNESDestroy() 23 */ 24 #undef __FUNCT__ 25 #define __FUNCT__ "SNESDestroy_NRichardson" 26 PetscErrorCode SNESDestroy_NRichardson(SNES snes) 27 { 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 ierr = SNESReset_NRichardson(snes);CHKERRQ(ierr); 32 ierr = PetscFree(snes->data);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 36 /* 37 SNESSetUp_NRichardson - Sets up the internal data structures for the later use 38 of the SNESNRICHARDSON nonlinear solver. 39 40 Input Parameters: 41 + snes - the SNES context 42 - x - the solution vector 43 44 Application Interface Routine: SNESSetUp() 45 */ 46 #undef __FUNCT__ 47 #define __FUNCT__ "SNESSetUp_NRichardson" 48 PetscErrorCode SNESSetUp_NRichardson(SNES snes) 49 { 50 const char *optionsprefix; 51 SNES_NRichardson *neP = (SNES_NRichardson *)snes->data; 52 PetscErrorCode ierr; 53 54 PetscFunctionBegin; 55 56 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 57 ierr = PetscLineSearchCreate(((PetscObject)snes)->comm, &neP->linesearch);CHKERRQ(ierr); 58 ierr = PetscLineSearchSetSNES(neP->linesearch, snes);CHKERRQ(ierr); 59 ierr = PetscLineSearchSetType(neP->linesearch, PETSCLINESEARCHL2);CHKERRQ(ierr); 60 ierr = PetscLineSearchAppendOptionsPrefix(neP->linesearch, optionsprefix);CHKERRQ(ierr); 61 ierr = PetscLineSearchSetFromOptions(neP->linesearch);CHKERRQ(ierr); 62 63 PetscFunctionReturn(0); 64 } 65 66 /* 67 SNESSetFromOptions_NRichardson - Sets various parameters for the SNESLS method. 68 69 Input Parameter: 70 . snes - the SNES context 71 72 Application Interface Routine: SNESSetFromOptions() 73 */ 74 #undef __FUNCT__ 75 #define __FUNCT__ "SNESSetFromOptions_NRichardson" 76 static PetscErrorCode SNESSetFromOptions_NRichardson(SNES snes) 77 { 78 PetscErrorCode ierr; 79 PetscFunctionBegin; 80 ierr = PetscOptionsHead("SNES Richardson options");CHKERRQ(ierr); 81 ierr = PetscOptionsTail();CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 /* 86 SNESView_NRichardson - Prints info from the SNESRichardson data structure. 87 88 Input Parameters: 89 + SNES - the SNES context 90 - viewer - visualization context 91 92 Application Interface Routine: SNESView() 93 */ 94 #undef __FUNCT__ 95 #define __FUNCT__ "SNESView_NRichardson" 96 static PetscErrorCode SNESView_NRichardson(SNES snes, PetscViewer viewer) 97 { 98 PetscBool iascii; 99 PetscErrorCode ierr; 100 SNES_NRichardson *neP = (SNES_NRichardson *)snes->data; 101 102 PetscFunctionBegin; 103 ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 104 if (iascii) { 105 ierr = PetscViewerASCIIPrintf(viewer," Line search variant: %s\n", ((PetscObject)neP->linesearch)->type_name);CHKERRQ(ierr); 106 } 107 PetscFunctionReturn(0); 108 } 109 110 /* 111 SNESSolve_NRichardson - Solves a nonlinear system with the Richardson method. 112 113 Input Parameters: 114 . snes - the SNES context 115 116 Output Parameter: 117 . outits - number of iterations until termination 118 119 Application Interface Routine: SNESSolve() 120 */ 121 #undef __FUNCT__ 122 #define __FUNCT__ "SNESSolve_NRichardson" 123 PetscErrorCode SNESSolve_NRichardson(SNES snes) 124 { 125 SNES_NRichardson *neP = (SNES_NRichardson *)snes->data; 126 Vec X, Y, F; 127 PetscReal xnorm, fnorm, ynorm; 128 PetscInt maxits, i; 129 PetscErrorCode ierr; 130 PetscBool lsSuccess; 131 SNESConvergedReason reason; 132 133 PetscFunctionBegin; 134 snes->reason = SNES_CONVERGED_ITERATING; 135 136 maxits = snes->max_its; /* maximum number of iterations */ 137 X = snes->vec_sol; /* X^n */ 138 Y = snes->vec_sol_update; /* \tilde X */ 139 F = snes->vec_func; /* residual vector */ 140 141 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 142 snes->iter = 0; 143 snes->norm = 0.; 144 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 145 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 146 if (snes->domainerror) { 147 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 148 PetscFunctionReturn(0); 149 } 150 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 151 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 152 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 153 snes->norm = fnorm; 154 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 155 SNESLogConvHistory(snes,fnorm,0); 156 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 157 158 /* set parameter for default relative tolerance convergence test */ 159 snes->ttol = fnorm*snes->rtol; 160 /* test convergence */ 161 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 162 if (snes->reason) PetscFunctionReturn(0); 163 164 for(i = 0; i < maxits; i++) { 165 lsSuccess = PETSC_TRUE; 166 /* Call general purpose update function */ 167 if (snes->ops->update) { 168 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 169 } 170 else if (snes->pc) { 171 ierr = VecCopy(X,Y);CHKERRQ(ierr); 172 ierr = SNESSolve(snes->pc, snes->vec_rhs, Y);CHKERRQ(ierr); 173 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 174 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 175 snes->reason = SNES_DIVERGED_INNER; 176 PetscFunctionReturn(0); 177 } 178 ierr = VecAYPX(Y,-1.0,X);CHKERRQ(ierr); 179 } else { 180 ierr = VecCopy(F,Y);CHKERRQ(ierr); 181 } 182 ierr = PetscLineSearchApply(neP->linesearch, X, F, &fnorm, Y);CHKERRQ(ierr); 183 ierr = PetscLineSearchGetNorms(neP->linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 184 ierr = PetscLineSearchGetSuccess(neP->linesearch, &lsSuccess);CHKERRQ(ierr); 185 if (!lsSuccess) { 186 if (++snes->numFailures >= snes->maxFailures) { 187 snes->reason = SNES_DIVERGED_LINE_SEARCH; 188 break; 189 } 190 } 191 if (snes->nfuncs >= snes->max_funcs) { 192 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 193 break; 194 } 195 if (snes->domainerror) { 196 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 197 PetscFunctionReturn(0); 198 } 199 /* Monitor convergence */ 200 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 201 snes->iter = i+1; 202 snes->norm = fnorm; 203 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 204 SNESLogConvHistory(snes,snes->norm,0); 205 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 206 /* Test for convergence */ 207 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 208 if (snes->reason) break; 209 } 210 if (i == maxits) { 211 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);CHKERRQ(ierr); 212 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 213 } 214 PetscFunctionReturn(0); 215 } 216 217 /*MC 218 SNESNRICHARDSON - Richardson nonlinear solver that uses successive substitutions, also sometimes known as Picard iteration. 219 220 Level: beginner 221 222 Options Database: 223 . -snes_ls <basic,basicnormnorms,quadratic,critical> Line search type. 224 225 Notes: If no inner nonlinear preconditioner is provided then solves F(x) - b = 0 using x^{n+1} = x^{n} - lambda 226 (F(x^n) - b) where lambda is obtained either SNESLineSearchSetDamping(), -snes_damping or a line search. If 227 an inner nonlinear preconditioner is provided (either with -npc_snes_type or SNESSetPC()) then the inner 228 solver is called an initial solution x^n and the nonlinear Richardson uses x^{n+1} = x^{n} + lambda d^{n} 229 where d^{n} = \hat{x}^{n} - x^{n} where \hat{x}^{n} is the solution returned from the inner solver. 230 231 This uses no derivative information thus will be much slower then Newton's method obtained with -snes_type ls 232 233 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESTR, SNESNGMRES, SNESQN, SNESNCG 234 M*/ 235 EXTERN_C_BEGIN 236 #undef __FUNCT__ 237 #define __FUNCT__ "SNESCreate_NRichardson" 238 PetscErrorCode SNESCreate_NRichardson(SNES snes) 239 { 240 PetscErrorCode ierr; 241 SNES_NRichardson *neP; 242 PetscFunctionBegin; 243 snes->ops->destroy = SNESDestroy_NRichardson; 244 snes->ops->setup = SNESSetUp_NRichardson; 245 snes->ops->setfromoptions = SNESSetFromOptions_NRichardson; 246 snes->ops->view = SNESView_NRichardson; 247 snes->ops->solve = SNESSolve_NRichardson; 248 snes->ops->reset = SNESReset_NRichardson; 249 250 snes->usesksp = PETSC_FALSE; 251 snes->usespc = PETSC_TRUE; 252 253 ierr = PetscNewLog(snes, SNES_NRichardson, &neP);CHKERRQ(ierr); 254 snes->data = (void*) neP; 255 256 snes->max_funcs = 30000; 257 snes->max_its = 10000; 258 259 PetscFunctionReturn(0); 260 } 261 EXTERN_C_END 262