1 #include <../src/snes/impls/gs/gsimpl.h> /*I "petscsnes.h" I*/ 2 3 /*@ 4 SNESNGSSetTolerances - Sets various parameters used in convergence tests. 5 6 Logically Collective on SNES 7 8 Input Parameters: 9 + snes - the SNES context 10 . abstol - absolute convergence tolerance 11 . rtol - relative convergence tolerance 12 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 13 - maxit - maximum number of iterations 14 15 Options Database Keys: 16 + -snes_ngs_atol <abstol> - Sets abstol 17 . -snes_ngs_rtol <rtol> - Sets rtol 18 . -snes_ngs_stol <stol> - Sets stol 19 - -snes_max_it <maxit> - Sets maxit 20 21 Level: intermediate 22 23 .seealso: SNESSetTrustRegionTolerance() 24 @*/ 25 PetscErrorCode SNESNGSSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit) 26 { 27 SNES_NGS *gs = (SNES_NGS*)snes->data; 28 29 PetscFunctionBegin; 30 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 31 32 if (abstol != PETSC_DEFAULT) { 33 if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol); 34 gs->abstol = abstol; 35 } 36 if (rtol != PETSC_DEFAULT) { 37 if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %g must be non-negative and less than 1.0",(double)rtol); 38 gs->rtol = rtol; 39 } 40 if (stol != PETSC_DEFAULT) { 41 if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %g must be non-negative",(double)stol); 42 gs->stol = stol; 43 } 44 if (maxit != PETSC_DEFAULT) { 45 if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 46 gs->max_its = maxit; 47 } 48 PetscFunctionReturn(0); 49 } 50 51 /*@ 52 SNESNGSGetTolerances - Gets various parameters used in convergence tests. 53 54 Not Collective 55 56 Input Parameters: 57 + snes - the SNES context 58 . atol - absolute convergence tolerance 59 . rtol - relative convergence tolerance 60 . stol - convergence tolerance in terms of the norm 61 of the change in the solution between steps 62 - maxit - maximum number of iterations 63 64 Notes: 65 The user can specify NULL for any parameter that is not needed. 66 67 Level: intermediate 68 69 .seealso: SNESSetTolerances() 70 @*/ 71 PetscErrorCode SNESNGSGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit) 72 { 73 SNES_NGS *gs = (SNES_NGS*)snes->data; 74 75 PetscFunctionBegin; 76 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 77 if (atol) *atol = gs->abstol; 78 if (rtol) *rtol = gs->rtol; 79 if (stol) *stol = gs->stol; 80 if (maxit) *maxit = gs->max_its; 81 PetscFunctionReturn(0); 82 } 83 84 /*@ 85 SNESNGSSetSweeps - Sets the number of sweeps of GS to use. 86 87 Input Parameters: 88 + snes - the SNES context 89 - sweeps - the number of sweeps of GS to perform. 90 91 Level: intermediate 92 93 .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetNPC(), SNESNGSGetSweeps() 94 @*/ 95 96 PetscErrorCode SNESNGSSetSweeps(SNES snes, PetscInt sweeps) 97 { 98 SNES_NGS *gs = (SNES_NGS*)snes->data; 99 100 PetscFunctionBegin; 101 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 102 gs->sweeps = sweeps; 103 PetscFunctionReturn(0); 104 } 105 106 /*@ 107 SNESNGSGetSweeps - Gets the number of sweeps GS will use. 108 109 Input Parameters: 110 . snes - the SNES context 111 112 Output Parameters: 113 . sweeps - the number of sweeps of GS to perform. 114 115 Level: intermediate 116 117 .seealso: SNESSetNGS(), SNESGetNGS(), SNESSetNPC(), SNESNGSSetSweeps() 118 @*/ 119 PetscErrorCode SNESNGSGetSweeps(SNES snes, PetscInt * sweeps) 120 { 121 SNES_NGS *gs = (SNES_NGS*)snes->data; 122 123 PetscFunctionBegin; 124 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 125 *sweeps = gs->sweeps; 126 PetscFunctionReturn(0); 127 } 128 129 PetscErrorCode SNESReset_NGS(SNES snes) 130 { 131 SNES_NGS *gs = (SNES_NGS*)snes->data; 132 PetscErrorCode ierr; 133 134 PetscFunctionBegin; 135 ierr = ISColoringDestroy(&gs->coloring);CHKERRQ(ierr); 136 PetscFunctionReturn(0); 137 } 138 139 PetscErrorCode SNESDestroy_NGS(SNES snes) 140 { 141 PetscErrorCode ierr; 142 143 PetscFunctionBegin; 144 ierr = SNESReset_NGS(snes);CHKERRQ(ierr); 145 ierr = PetscFree(snes->data);CHKERRQ(ierr); 146 PetscFunctionReturn(0); 147 } 148 149 PetscErrorCode SNESSetUp_NGS(SNES snes) 150 { 151 PetscErrorCode ierr; 152 PetscErrorCode (*f)(SNES,Vec,Vec,void*); 153 154 PetscFunctionBegin; 155 ierr = SNESGetNGS(snes,&f,NULL);CHKERRQ(ierr); 156 if (!f) { 157 ierr = SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);CHKERRQ(ierr); 158 } 159 PetscFunctionReturn(0); 160 } 161 162 PetscErrorCode SNESSetFromOptions_NGS(PetscOptionItems *PetscOptionsObject,SNES snes) 163 { 164 SNES_NGS *gs = (SNES_NGS*)snes->data; 165 PetscErrorCode ierr; 166 PetscInt sweeps,max_its=PETSC_DEFAULT; 167 PetscReal rtol=PETSC_DEFAULT,atol=PETSC_DEFAULT,stol=PETSC_DEFAULT; 168 PetscBool flg,flg1,flg2,flg3; 169 170 PetscFunctionBegin; 171 ierr = PetscOptionsHead(PetscOptionsObject,"SNES GS options");CHKERRQ(ierr); 172 /* GS Options */ 173 ierr = PetscOptionsInt("-snes_ngs_sweeps","Number of sweeps of GS to apply","SNESComputeGS",gs->sweeps,&sweeps,&flg);CHKERRQ(ierr); 174 if (flg) { 175 ierr = SNESNGSSetSweeps(snes,sweeps);CHKERRQ(ierr); 176 } 177 ierr = PetscOptionsReal("-snes_ngs_atol","Absolute residual tolerance for GS iteration","SNESComputeGS",gs->abstol,&atol,&flg);CHKERRQ(ierr); 178 ierr = PetscOptionsReal("-snes_ngs_rtol","Relative residual tolerance for GS iteration","SNESComputeGS",gs->rtol,&rtol,&flg1);CHKERRQ(ierr); 179 ierr = PetscOptionsReal("-snes_ngs_stol","Absolute update tolerance for GS iteration","SNESComputeGS",gs->stol,&stol,&flg2);CHKERRQ(ierr); 180 ierr = PetscOptionsInt("-snes_ngs_max_it","Maximum number of sweeps of GS to apply","SNESComputeGS",gs->max_its,&max_its,&flg3);CHKERRQ(ierr); 181 if (flg || flg1 || flg2 || flg3) { 182 ierr = SNESNGSSetTolerances(snes,atol,rtol,stol,max_its);CHKERRQ(ierr); 183 } 184 flg = PETSC_FALSE; 185 ierr = PetscOptionsBool("-snes_ngs_secant","Use finite difference secant approximation with coloring","",flg,&flg,NULL);CHKERRQ(ierr); 186 if (flg) { 187 ierr = SNESSetNGS(snes,SNESComputeNGSDefaultSecant,NULL);CHKERRQ(ierr); 188 ierr = PetscInfo(snes,"Setting default finite difference secant approximation with coloring\n");CHKERRQ(ierr); 189 } 190 ierr = PetscOptionsReal("-snes_ngs_secant_h","Differencing parameter for secant search","",gs->h,&gs->h,NULL);CHKERRQ(ierr); 191 ierr = PetscOptionsBool("-snes_ngs_secant_mat_coloring","Use the graph coloring of the Jacobian for the secant GS","",gs->secant_mat,&gs->secant_mat,&flg);CHKERRQ(ierr); 192 193 ierr = PetscOptionsTail();CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196 197 PetscErrorCode SNESView_NGS(SNES snes, PetscViewer viewer) 198 { 199 PetscErrorCode ierr; 200 PetscErrorCode (*f)(SNES,Vec,Vec,void*); 201 SNES_NGS *gs = (SNES_NGS*)snes->data; 202 PetscBool iascii; 203 204 PetscFunctionBegin; 205 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 206 if (iascii) { 207 ierr = DMSNESGetNGS(snes->dm,&f,NULL);CHKERRQ(ierr); 208 if (f == SNESComputeNGSDefaultSecant) { 209 ierr = PetscViewerASCIIPrintf(viewer," Use finite difference secant approximation with coloring with h = %g \n",(double)gs->h);CHKERRQ(ierr); 210 } 211 } 212 PetscFunctionReturn(0); 213 } 214 215 PetscErrorCode SNESSolve_NGS(SNES snes) 216 { 217 Vec F; 218 Vec X; 219 Vec B; 220 PetscInt i; 221 PetscReal fnorm; 222 PetscErrorCode ierr; 223 SNESNormSchedule normschedule; 224 225 PetscFunctionBegin; 226 227 if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name); 228 229 ierr = PetscCitationsRegister(SNESCitation,&SNEScite);CHKERRQ(ierr); 230 X = snes->vec_sol; 231 F = snes->vec_func; 232 B = snes->vec_rhs; 233 234 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 235 snes->iter = 0; 236 snes->norm = 0.; 237 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 238 snes->reason = SNES_CONVERGED_ITERATING; 239 240 ierr = SNESGetNormSchedule(snes, &normschedule);CHKERRQ(ierr); 241 if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) { 242 /* compute the initial function and preconditioned update delX */ 243 if (!snes->vec_func_init_set) { 244 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 245 } else snes->vec_func_init_set = PETSC_FALSE; 246 247 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 248 SNESCheckFunctionNorm(snes,fnorm); 249 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 250 snes->iter = 0; 251 snes->norm = fnorm; 252 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 253 ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr); 254 ierr = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr); 255 256 /* test convergence */ 257 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 258 if (snes->reason) PetscFunctionReturn(0); 259 } else { 260 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 261 ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr); 262 } 263 264 /* Call general purpose update function */ 265 if (snes->ops->update) { 266 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 267 } 268 269 for (i = 0; i < snes->max_its; i++) { 270 ierr = SNESComputeNGS(snes, B, X);CHKERRQ(ierr); 271 /* only compute norms if requested or about to exit due to maximum iterations */ 272 if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) { 273 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 274 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 275 SNESCheckFunctionNorm(snes,fnorm); 276 /* Monitor convergence */ 277 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 278 snes->iter = i+1; 279 snes->norm = fnorm; 280 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 281 ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr); 282 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 283 } 284 /* Test for convergence */ 285 if (normschedule == SNES_NORM_ALWAYS) {ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);} 286 if (snes->reason) PetscFunctionReturn(0); 287 /* Call general purpose update function */ 288 if (snes->ops->update) { 289 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 290 } 291 } 292 if (normschedule == SNES_NORM_ALWAYS) { 293 if (i == snes->max_its) { 294 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);CHKERRQ(ierr); 295 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 296 } 297 } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* GS is meant to be used as a preconditioner */ 298 PetscFunctionReturn(0); 299 } 300 301 /*MC 302 SNESNGS - Either calls the user-provided solution routine provided with SNESSetNGS() or does a finite difference secant approximation 303 using coloring. 304 305 Level: advanced 306 307 Options Database: 308 + -snes_ngs_sweeps <n> - Number of sweeps of GS to apply 309 . -snes_ngs_atol <atol> - Absolute residual tolerance for GS iteration 310 . -snes_ngs_rtol <rtol> - Relative residual tolerance for GS iteration 311 . -snes_ngs_stol <stol> - Absolute update tolerance for GS iteration 312 . -snes_ngs_max_it <maxit> - Maximum number of sweeps of GS to apply 313 . -snes_ngs_secant - Use pointwise secant local Jacobian approximation with coloring instead of user provided Gauss-Seidel routine, this is 314 used by default if no user provided Gauss-Seidel routine is available. Requires either that a DM that can compute a coloring 315 is available or a Jacobian sparse matrix is provided (from which to get the coloring). 316 . -snes_ngs_secant_h <h> - Differencing parameter for secant approximation 317 . -snes_ngs_secant_mat_coloring - Use the graph coloring of the Jacobian for the secant GS even if a DM is available. 318 - -snes_norm_schedule <none, always, initialonly, finalonly, initialfinalonly> - how often the residual norms are computed 319 320 Notes: 321 the Gauss-Seidel smoother is inherited through composition. If a solver has been created with SNESGetNPC(), it will have 322 its parent's Gauss-Seidel routine associated with it. 323 324 By default this routine computes the solution norm at each iteration, this can be time consuming, you can turn this off with SNESSetNormSchedule() or -snes_norm_schedule 325 References: 326 . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers", 327 SIAM Review, 57(4), 2015 328 329 .seealso: SNESCreate(), SNES, SNESSetType(), SNESSetNGS(), SNESType (for list of available types), SNESNGSSetSweeps(), SNESNGSSetTolerances(), 330 SNESSetNormSchedule() 331 M*/ 332 333 PETSC_EXTERN PetscErrorCode SNESCreate_NGS(SNES snes) 334 { 335 SNES_NGS *gs; 336 PetscErrorCode ierr; 337 338 PetscFunctionBegin; 339 snes->ops->destroy = SNESDestroy_NGS; 340 snes->ops->setup = SNESSetUp_NGS; 341 snes->ops->setfromoptions = SNESSetFromOptions_NGS; 342 snes->ops->view = SNESView_NGS; 343 snes->ops->solve = SNESSolve_NGS; 344 snes->ops->reset = SNESReset_NGS; 345 346 snes->usesksp = PETSC_FALSE; 347 snes->usesnpc = PETSC_FALSE; 348 349 snes->alwayscomputesfinalresidual = PETSC_FALSE; 350 351 if (!snes->tolerancesset) { 352 snes->max_its = 10000; 353 snes->max_funcs = 10000; 354 } 355 356 ierr = PetscNewLog(snes,&gs);CHKERRQ(ierr); 357 358 gs->sweeps = 1; 359 gs->rtol = 1e-5; 360 gs->abstol = PETSC_MACHINE_EPSILON; 361 gs->stol = 1000*PETSC_MACHINE_EPSILON; 362 gs->max_its = 50; 363 gs->h = PETSC_SQRT_MACHINE_EPSILON; 364 365 snes->data = (void*) gs; 366 PetscFunctionReturn(0); 367 } 368