1 #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/ 2 3 typedef struct { 4 PetscInt sweeps; 5 } SNES_GS; 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "SNESReset_GS" 9 PetscErrorCode SNESReset_GS(SNES snes) 10 { 11 PetscFunctionBegin; 12 PetscFunctionReturn(0); 13 } 14 15 #undef __FUNCT__ 16 #define __FUNCT__ "SNESDestroy_GS" 17 PetscErrorCode SNESDestroy_GS(SNES snes) 18 { 19 PetscErrorCode ierr; 20 21 PetscFunctionBegin; 22 ierr = SNESReset_GS(snes);CHKERRQ(ierr); 23 ierr = PetscFree(snes->data); 24 PetscFunctionReturn(0); 25 } 26 27 #undef __FUNCT__ 28 #define __FUNCT__ "SNESSetUp_GS" 29 PetscErrorCode SNESSetUp_GS(SNES snes) 30 { 31 PetscFunctionBegin; 32 PetscFunctionReturn(0); 33 } 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "SNESSetFromOptions_GS" 37 PetscErrorCode SNESSetFromOptions_GS(SNES snes) 38 { 39 PetscErrorCode ierr; 40 41 PetscFunctionBegin; 42 ierr = PetscOptionsHead("SNES GS options");CHKERRQ(ierr); 43 ierr = PetscOptionsTail();CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "SNESView_GS" 49 PetscErrorCode SNESView_GS(SNES snes, PetscViewer viewer) 50 { 51 PetscFunctionBegin; 52 PetscFunctionReturn(0); 53 } 54 55 #undef __FUNCT__ 56 #define __FUNCT__ "SNESSolve_GS" 57 PetscErrorCode SNESSolve_GS(SNES snes) 58 { 59 Vec F; 60 Vec X; 61 Vec B; 62 PetscInt i; 63 PetscReal fnorm; 64 PetscErrorCode ierr; 65 SNESNormType normtype; 66 67 PetscFunctionBegin; 68 69 X = snes->vec_sol; 70 F = snes->vec_func; 71 B = snes->vec_rhs; 72 73 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 74 snes->iter = 0; 75 snes->norm = 0.; 76 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 77 snes->reason = SNES_CONVERGED_ITERATING; 78 79 ierr = SNESGetNormType(snes, &normtype);CHKERRQ(ierr); 80 if (normtype == SNES_NORM_FUNCTION || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) { 81 /* compute the initial function and preconditioned update delX */ 82 if (!snes->vec_func_init_set) { 83 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 84 if (snes->domainerror) { 85 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 86 PetscFunctionReturn(0); 87 } 88 } else { 89 snes->vec_func_init_set = PETSC_FALSE; 90 } 91 92 /* convergence test */ 93 if (!snes->norm_init_set) { 94 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 95 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 96 } else { 97 fnorm = snes->norm_init; 98 snes->norm_init_set = PETSC_FALSE; 99 } 100 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 101 snes->iter = 0; 102 snes->norm = fnorm; 103 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 104 SNESLogConvHistory(snes,snes->norm,0); 105 ierr = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr); 106 107 /* set parameter for default relative tolerance convergence test */ 108 snes->ttol = fnorm*snes->rtol; 109 110 /* test convergence */ 111 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 112 if (snes->reason) PetscFunctionReturn(0); 113 } else { 114 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 115 SNESLogConvHistory(snes,snes->norm,0); 116 ierr = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr); 117 } 118 119 120 /* Call general purpose update function */ 121 if (snes->ops->update) { 122 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 123 } 124 125 for (i = 0; i < snes->max_its; i++) { 126 ierr = SNESComputeGS(snes, B, X);CHKERRQ(ierr); 127 /* only compute norms if requested or about to exit due to maximum iterations */ 128 if (normtype == SNES_NORM_FUNCTION || ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY))) { 129 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 130 if (snes->domainerror) { 131 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 132 PetscFunctionReturn(0); 133 } 134 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 135 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 136 } 137 /* Monitor convergence */ 138 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 139 snes->iter = i+1; 140 snes->norm = fnorm; 141 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 142 SNESLogConvHistory(snes,snes->norm,0); 143 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 144 /* Test for convergence */ 145 if (normtype == SNES_NORM_FUNCTION) ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 146 if (snes->reason) PetscFunctionReturn(0); 147 /* Call general purpose update function */ 148 if (snes->ops->update) { 149 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 150 } 151 } 152 if (normtype == SNES_NORM_FUNCTION) { 153 if (i == snes->max_its) { 154 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);CHKERRQ(ierr); 155 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 156 } 157 } else { 158 if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* GS is meant to be used as a preconditioner */ 159 } 160 PetscFunctionReturn(0); 161 } 162 163 /*MC 164 SNESGS - Just calls the user-provided solution routine provided with SNESSetGS() 165 166 Level: advanced 167 168 Notes: 169 the Gauss-Seidel smoother is inherited through composition. If a solver has been created with SNESGetPC(), it will have 170 its parent's Gauss-Seidel routine associated with it. 171 172 .seealso: SNESCreate(), SNES, SNESSetType(), SNESSetGS(), SNESType (for list of available types) 173 M*/ 174 175 EXTERN_C_BEGIN 176 #undef __FUNCT__ 177 #define __FUNCT__ "SNESCreate_GS" 178 PetscErrorCode SNESCreate_GS(SNES snes) 179 { 180 SNES_GS *gs; 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 snes->ops->destroy = SNESDestroy_GS; 185 snes->ops->setup = SNESSetUp_GS; 186 snes->ops->setfromoptions = SNESSetFromOptions_GS; 187 snes->ops->view = SNESView_GS; 188 snes->ops->solve = SNESSolve_GS; 189 snes->ops->reset = SNESReset_GS; 190 191 snes->usesksp = PETSC_FALSE; 192 snes->usespc = PETSC_FALSE; 193 194 if (!snes->tolerancesset) { 195 snes->max_its = 10000; 196 snes->max_funcs = 10000; 197 } 198 199 ierr = PetscNewLog(snes, SNES_GS, &gs);CHKERRQ(ierr); 200 snes->data = (void*) gs; 201 PetscFunctionReturn(0); 202 } 203 EXTERN_C_END 204