1 /* Defines the basic SNES object */ 2 #include <private/snesimpl.h> 3 4 /* Private structure for the Anderson mixing method aka nonlinear Krylov */ 5 typedef struct { 6 Vec *v, *w; 7 PetscReal *f2; /* 2-norms of function (residual) at each stage */ 8 PetscInt msize; /* maximum size of space */ 9 PetscInt csize; /* current size of space */ 10 } SNES_NGMRES; 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "SNESReset_NGMRES" 14 PetscErrorCode SNESReset_NGMRES(SNES snes) 15 { 16 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 17 PetscErrorCode ierr; 18 19 PetscFunctionBegin; 20 ierr = VecDestroyVecs(ngmres->msize, &ngmres->v);CHKERRQ(ierr); 21 ierr = VecDestroyVecs(ngmres->msize, &ngmres->w);CHKERRQ(ierr); 22 PetscFunctionReturn(0); 23 } 24 25 #undef __FUNCT__ 26 #define __FUNCT__ "SNESDestroy_NGMRES" 27 PetscErrorCode SNESDestroy_NGMRES(SNES snes) 28 { 29 PetscErrorCode ierr; 30 31 PetscFunctionBegin; 32 ierr = SNESReset_NGMRES(snes);CHKERRQ(ierr); 33 if (snes->work) {ierr = VecDestroyVecs(snes->nwork, &snes->work);CHKERRQ(ierr);} 34 PetscFunctionReturn(0); 35 } 36 37 #undef __FUNCT__ 38 #define __FUNCT__ "SNESSetUp_NGMRES" 39 PetscErrorCode SNESSetUp_NGMRES(SNES snes) 40 { 41 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 42 PetscErrorCode ierr; 43 44 PetscFunctionBegin; 45 #if 0 46 if (snes->pc_side != PC_LEFT) {SETERRQ(((PetscObject) snes)->comm, PETSC_ERR_SUP, "Only left preconditioning allowed for SNESNGMRES");} 47 #endif 48 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->v);CHKERRQ(ierr); 49 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->w);CHKERRQ(ierr); 50 ierr = SNESDefaultGetWork(snes, 1);CHKERRQ(ierr); 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNCT__ 55 #define __FUNCT__ "SNESSetFromOptions_NGMRES" 56 PetscErrorCode SNESSetFromOptions_NGMRES(SNES snes) 57 { 58 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 59 PetscErrorCode ierr; 60 61 PetscFunctionBegin; 62 ierr = PetscOptionsHead("SNES NGMRES options");CHKERRQ(ierr); 63 ierr = PetscOptionsInt("-snes_gmres_restart", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, PETSC_NULL);CHKERRQ(ierr); 64 ierr = PetscOptionsTail();CHKERRQ(ierr); 65 PetscFunctionReturn(0); 66 } 67 68 #undef __FUNCT__ 69 #define __FUNCT__ "SNESView_NGMRES" 70 PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer) 71 { 72 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 73 PetscBool iascii; 74 PetscErrorCode ierr; 75 76 PetscFunctionBegin; 77 ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 78 if (iascii) { 79 ierr = PetscViewerASCIIPrintf(viewer, " Size of space %d\n", ngmres->msize);CHKERRQ(ierr); 80 } else { 81 SETERRQ1(((PetscObject) snes)->comm, PETSC_ERR_SUP, "Viewer type %s not supported for SNESNGMRES", ((PetscObject) viewer)->type_name); 82 } 83 PetscFunctionReturn(0); 84 } 85 86 #undef __FUNCT__ 87 #define __FUNCT__ "SNESSolve_NGMRES" 88 PetscErrorCode SNESSolve_NGMRES(SNES snes) 89 { 90 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 91 SNES pc; 92 Vec X, Y, F, r, rOld, *V = ngmres->v, *W = ngmres->w; 93 PetscScalar wdot; 94 PetscReal fnorm; 95 PetscInt i, j, k; 96 PetscErrorCode ierr; 97 98 PetscFunctionBegin; 99 snes->reason = SNES_CONVERGED_ITERATING; 100 X = snes->vec_sol; 101 Y = snes->vec_sol_update; 102 F = snes->vec_func; 103 rOld = snes->work[0]; 104 105 ierr = SNESGetPC(snes, &pc);CHKERRQ(ierr); 106 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 107 snes->iter = 0; 108 snes->norm = 0.; 109 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 110 ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr); 111 if (snes->domainerror) { 112 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 113 PetscFunctionReturn(0); 114 } 115 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); 116 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in norm"); 117 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 118 snes->norm = fnorm; 119 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 120 SNESLogConvHistory(snes, fnorm, 0); 121 ierr = SNESMonitor(snes, 0, fnorm);CHKERRQ(ierr); 122 123 /* set parameter for default relative tolerance convergence test */ 124 snes->ttol = fnorm*snes->rtol; 125 /* test convergence */ 126 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 127 if (snes->reason) PetscFunctionReturn(0); 128 129 #ifndef OLD_BARRY_CODE 130 ierr = VecCopy(F, rOld);CHKERRQ(ierr); 131 #else 132 /* Barry: What the heck is this part doing? */ 133 /* determine optimal scale factor -- slow code */ 134 ierr = VecDuplicate(P, &y);CHKERRQ(ierr); 135 ierr = VecDuplicate(P, &w);CHKERRQ(ierr); 136 ierr = MatMult(Amat, Pold, y);CHKERRQ(ierr); 137 /*ierr = KSP_PCApplyBAorAB(ksp,Pold,y,w);CHKERRQ(ierr); */ /* y = BAp */ 138 ierr = VecDotNorm2(Pold, y, &rdot, &abr);CHKERRQ(ierr); /* rdot = (p)^T(BAp); abr = (BAp)^T (BAp) */ 139 ierr = VecDestroy(&y);CHKERRQ(ierr); 140 ierr = VecDestroy(&w);CHKERRQ(ierr); 141 A0 = rdot/abr; 142 ierr = VecAXPY(X,A0,Pold);CHKERRQ(ierr); /* x <- x + scale p */ 143 #endif 144 145 /* Loop over batches of directions */ 146 /* Code from Barry I do not understand to solve the least-squares problem, Time to try again */ 147 for(k = 0; k < snes->max_its; k += ngmres->msize) { 148 /* Loop over updates for this batch */ 149 /* TODO: Incorporate the variant which use the analytic Jacobian */ 150 /* TODO: Incorporate selection criteria for u^A from paper (need to store residual and update norms) */ 151 /* TODO: Incorporate criteria for restarting from paper */ 152 for(i = 0; i < ngmres->msize && k+i < snes->max_its; ++i) { 153 /* Get a new approxiate solution u^k */ 154 ierr = SNESSolve(pc, PETSC_NULL, X);CHKERRQ(ierr); 155 /* Solve least squares problem using last msize iterates */ 156 ierr = SNESGetFunction(pc, &r, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); /* r = F(u^M) */ 157 for(j = 0; j < i; ++j) { /* r = \prod_j (I + v_j w^T_j) r */ 158 ierr = VecDot(W[j], r, &wdot);CHKERRQ(ierr); 159 ierr = VecAXPY(r, wdot, V[j]);CHKERRQ(ierr); 160 } 161 ierr = VecCopy(rOld, W[i]);CHKERRQ(ierr); /* w_i = r_{old} */ 162 ierr = VecAXPY(rOld, -1.0, r);CHKERRQ(ierr); /* v_i = r */ 163 ierr = VecDot(W[i], rOld, &wdot);CHKERRQ(ierr); /* ------------------- */ 164 ierr = VecCopy(r, V[i]);CHKERRQ(ierr); /* w^T_i (r_{old} - r) */ 165 ierr = VecScale(V[i], 1.0/wdot);CHKERRQ(ierr); 166 ierr = VecDot(W[i], r, &wdot);CHKERRQ(ierr); /* r = (I + v_i w^T_i) r */ 167 ierr = VecAXPY(r, wdot, V[i]);CHKERRQ(ierr); 168 ierr = VecCopy(r, rOld);CHKERRQ(ierr); /* r_{old} = r */ 169 170 ierr = VecAXPY(X, 1.0, r);CHKERRQ(ierr); /* u^A = u^M + r */ 171 /* Monitor and log history */ 172 ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr); /* F(u^k) */ 173 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); 174 SNESLogConvHistory(snes, fnorm, 0); 175 ierr = SNESMonitor(snes, i+k+1, fnorm);CHKERRQ(ierr); 176 /* Check convergence */ 177 ierr = (*snes->ops->converged)(snes, i+k+1, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);CHKERRQ(ierr); 178 if (snes->reason) PetscFunctionReturn(0); 179 /* Select update between u^M and u^A */ 180 /* Decide whether to restart */ 181 } 182 } 183 snes->reason = SNES_DIVERGED_MAX_IT; 184 PetscFunctionReturn(0); 185 } 186 187 /*MC 188 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 189 190 Level: beginner 191 192 Notes: Supports only left preconditioning 193 194 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 195 SIAM Journal on Scientific Computing, 21(5), 2000. 196 197 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 198 M*/ 199 EXTERN_C_BEGIN 200 #undef __FUNCT__ 201 #define __FUNCT__ "SNESCreate_NGMRES" 202 PetscErrorCode SNESCreate_NGMRES(SNES snes) 203 { 204 SNES_NGMRES *ngmres; 205 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 snes->ops->destroy = SNESDestroy_NGMRES; 209 snes->ops->setup = SNESSetUp_NGMRES; 210 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 211 snes->ops->view = SNESView_NGMRES; 212 snes->ops->solve = SNESSolve_NGMRES; 213 snes->ops->reset = SNESReset_NGMRES; 214 215 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 216 snes->data = (void*) ngmres; 217 ngmres->msize = 30; 218 ngmres->csize = 0; 219 220 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 221 #if 0 222 if (ksp->pc_side != PC_LEFT) {ierr = PetscInfo(ksp,"WARNING! Setting PC_SIDE for NGMRES to left!\n");CHKERRQ(ierr);} 223 snes->pc_side = PC_LEFT; 224 #endif 225 PetscFunctionReturn(0); 226 } 227 EXTERN_C_END 228