Actual source code: anderson.c

      1: #include <../src/snes/impls/ngmres/snesngmres.h>

      3: static PetscErrorCode SNESSetFromOptions_Anderson(SNES snes, PetscOptionItems *PetscOptionsObject)
      4: {
      5:   SNES_NGMRES *ngmres  = (SNES_NGMRES *)snes->data;
      6:   PetscBool    monitor = PETSC_FALSE;

      8:   PetscFunctionBegin;
      9:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES NGMRES options");
     10:   PetscCall(PetscOptionsInt("-snes_anderson_m", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, NULL));
     11:   PetscCall(PetscOptionsReal("-snes_anderson_beta", "Mixing parameter", "SNES", ngmres->andersonBeta, &ngmres->andersonBeta, NULL));
     12:   PetscCall(PetscOptionsInt("-snes_anderson_restart", "Iterations before forced restart", "SNES", ngmres->restart_periodic, &ngmres->restart_periodic, NULL));
     13:   PetscCall(PetscOptionsInt("-snes_anderson_restart_it", "Tolerance iterations before restart", "SNES", ngmres->restart_it, &ngmres->restart_it, NULL));
     14:   PetscCall(PetscOptionsEnum("-snes_anderson_restart_type", "Restart type", "SNESNGMRESSetRestartType", SNESNGMRESRestartTypes, (PetscEnum)ngmres->restart_type, (PetscEnum *)&ngmres->restart_type, NULL));
     15:   PetscCall(PetscOptionsBool("-snes_anderson_monitor", "Monitor steps of Anderson Mixing", "SNES", ngmres->monitor ? PETSC_TRUE : PETSC_FALSE, &monitor, NULL));
     16:   if (monitor) ngmres->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
     17:   PetscOptionsHeadEnd();
     18:   PetscFunctionReturn(PETSC_SUCCESS);
     19: }

     21: static PetscErrorCode SNESSolve_Anderson(SNES snes)
     22: {
     23:   SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
     24:   /* present solution, residual, and preconditioned residual */
     25:   Vec X, F, B, D;
     26:   /* candidate linear combination answers */
     27:   Vec XA, FA, XM, FM;

     29:   /* coefficients and RHS to the minimization problem */
     30:   PetscReal           fnorm, fMnorm, fAnorm;
     31:   PetscReal           xnorm, ynorm;
     32:   PetscReal           dnorm, dminnorm = 0.0, fminnorm;
     33:   PetscInt            restart_count = 0;
     34:   PetscInt            k, k_restart, l, ivec;
     35:   PetscBool           selectRestart;
     36:   SNESConvergedReason reason;

     38:   PetscFunctionBegin;
     39:   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);

     41:   PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
     42:   /* variable initialization */
     43:   snes->reason = SNES_CONVERGED_ITERATING;
     44:   X            = snes->vec_sol;
     45:   F            = snes->vec_func;
     46:   B            = snes->vec_rhs;
     47:   XA           = snes->work[2];
     48:   FA           = snes->work[0];
     49:   D            = snes->work[1];

     51:   /* work for the line search */
     52:   XM = snes->work[3];
     53:   FM = snes->work[4];

     55:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
     56:   snes->iter = 0;
     57:   snes->norm = 0.;
     58:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

     60:   /* initialization */

     62:   /* r = F(x) */

     64:   if (snes->npc && snes->npcside == PC_LEFT) {
     65:     PetscCall(SNESApplyNPC(snes, X, NULL, F));
     66:     PetscCall(SNESGetConvergedReason(snes->npc, &reason));
     67:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
     68:       snes->reason = SNES_DIVERGED_INNER;
     69:       PetscFunctionReturn(PETSC_SUCCESS);
     70:     }
     71:     PetscCall(VecNorm(F, NORM_2, &fnorm));
     72:   } else {
     73:     if (!snes->vec_func_init_set) {
     74:       PetscCall(SNESComputeFunction(snes, X, F));
     75:     } else snes->vec_func_init_set = PETSC_FALSE;

     77:     PetscCall(VecNorm(F, NORM_2, &fnorm));
     78:     SNESCheckFunctionNorm(snes, fnorm);
     79:   }
     80:   fminnorm = fnorm;

     82:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
     83:   snes->norm = fnorm;
     84:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
     85:   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
     86:   PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
     87:   PetscCall(SNESMonitor(snes, 0, fnorm));
     88:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
     89:   PetscCall(SNESNGMRESUpdateSubspace_Private(snes, 0, 0, F, fnorm, X));

     91:   k_restart = 0;
     92:   l         = 0;
     93:   ivec      = 0;
     94:   for (k = 1; k < snes->max_its + 1; k++) {
     95:     /* Call general purpose update function */
     96:     PetscTryTypeMethod(snes, update, snes->iter);

     98:     /* select which vector of the stored subspace will be updated */
     99:     if (snes->npc && snes->npcside == PC_RIGHT) {
    100:       PetscCall(VecCopy(X, XM));
    101:       PetscCall(SNESSetInitialFunction(snes->npc, F));

    103:       PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, XM, B, 0));
    104:       PetscCall(SNESSolve(snes->npc, B, XM));
    105:       PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, XM, B, 0));

    107:       PetscCall(SNESGetConvergedReason(snes->npc, &reason));
    108:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
    109:         snes->reason = SNES_DIVERGED_INNER;
    110:         PetscFunctionReturn(PETSC_SUCCESS);
    111:       }
    112:       PetscCall(SNESGetNPCFunction(snes, FM, &fMnorm));
    113:       PetscCall(VecAXPBY(XM, 1.0 - ngmres->andersonBeta, ngmres->andersonBeta, X));
    114:     } else {
    115:       PetscCall(VecCopy(F, FM));
    116:       PetscCall(VecWAXPY(XM, -ngmres->andersonBeta, FM, X));
    117:       fMnorm = fnorm;
    118:     }

    120:     PetscCall(SNESNGMRESFormCombinedSolution_Private(snes, ivec, l, XM, FM, fMnorm, X, XA, FA));
    121:     ivec = k_restart % ngmres->msize;
    122:     if (ngmres->restart_type == SNES_NGMRES_RESTART_DIFFERENCE) {
    123:       PetscCall(SNESNGMRESNorms_Private(snes, l, X, F, XM, FM, XA, FA, D, &dnorm, &dminnorm, NULL, NULL, NULL, &xnorm, &fAnorm, &ynorm));
    124:       PetscCall(SNESNGMRESSelectRestart_Private(snes, l, snes->norm, fMnorm, fnorm, dnorm, fminnorm, dminnorm, &selectRestart));
    125:       /* if the restart conditions persist for more than restart_it iterations, restart. */
    126:       if (selectRestart) restart_count++;
    127:       else restart_count = 0;
    128:     } else if (ngmres->restart_type == SNES_NGMRES_RESTART_PERIODIC) {
    129:       PetscCall(SNESNGMRESNorms_Private(snes, l, X, F, XM, FM, XA, FA, D, NULL, NULL, NULL, NULL, NULL, &xnorm, &fAnorm, &ynorm));
    130:       if (k_restart > ngmres->restart_periodic) {
    131:         if (ngmres->monitor) PetscCall(PetscViewerASCIIPrintf(ngmres->monitor, "periodic restart after %" PetscInt_FMT " iterations\n", k_restart));
    132:         restart_count = ngmres->restart_it;
    133:       }
    134:     } else {
    135:       PetscCall(SNESNGMRESNorms_Private(snes, l, X, F, XM, FM, XA, FA, D, NULL, NULL, NULL, NULL, NULL, &xnorm, &fAnorm, &ynorm));
    136:     }
    137:     /* restart after restart conditions have persisted for a fixed number of iterations */
    138:     if (restart_count >= ngmres->restart_it) {
    139:       if (ngmres->monitor) PetscCall(PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %" PetscInt_FMT "\n", k_restart));
    140:       restart_count = 0;
    141:       k_restart     = 0;
    142:       l             = 0;
    143:       ivec          = 0;
    144:     } else {
    145:       if (l < ngmres->msize) l++;
    146:       k_restart++;
    147:       PetscCall(SNESNGMRESUpdateSubspace_Private(snes, ivec, l, FM, fnorm, XM));
    148:     }

    150:     fnorm = fAnorm;
    151:     if (fminnorm > fnorm) fminnorm = fnorm;

    153:     PetscCall(VecCopy(XA, X));
    154:     PetscCall(VecCopy(FA, F));

    156:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    157:     snes->iter  = k;
    158:     snes->norm  = fnorm;
    159:     snes->xnorm = xnorm;
    160:     snes->ynorm = ynorm;
    161:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
    162:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, snes->iter));
    163:     PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
    164:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
    165:     if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
    166:   }
    167:   snes->reason = SNES_DIVERGED_MAX_IT;
    168:   PetscFunctionReturn(PETSC_SUCCESS);
    169: }

    171: /*MC
    172:   SNESANDERSON - Implements the Anderson Mixing nonlinear solver {cite}`anderson1965`, {cite}`bruneknepleysmithtu15`

    174:    Level: beginner

    176:    Options Database Keys:
    177: +  -snes_anderson_m                - Number of stored previous solutions and residuals
    178: .  -snes_anderson_beta             - Anderson mixing parameter
    179: .  -snes_anderson_restart_type     - Type of restart (see `SNESNGMRES`)
    180: .  -snes_anderson_restart_it       - Number of iterations of restart conditions before restart
    181: .  -snes_anderson_restart          - Number of iterations before periodic restart
    182: -  -snes_anderson_monitor          - Prints relevant information about the Anderson mixing iteration

    184:    Notes:
    185:    The Anderson Mixing method combines m previous solutions into a minimum-residual solution by solving a small linearized
    186:    optimization problem at each iteration.

    188:    Very similar to the `SNESNGMRES` algorithm.

    190:    This algorithm ignores any Jacobian provided with `SNESSetJacobian()`

    192:    Only supports left non-linear preconditioning.

    194: .seealso: [](ch_snes), `SNESNGMRES`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESType`
    195: M*/

    197: PETSC_EXTERN PetscErrorCode SNESCreate_Anderson(SNES snes)
    198: {
    199:   SNES_NGMRES   *ngmres;
    200:   SNESLineSearch linesearch;

    202:   PetscFunctionBegin;
    203:   snes->ops->destroy        = SNESDestroy_NGMRES;
    204:   snes->ops->setup          = SNESSetUp_NGMRES;
    205:   snes->ops->setfromoptions = SNESSetFromOptions_Anderson;
    206:   snes->ops->view           = SNESView_NGMRES;
    207:   snes->ops->solve          = SNESSolve_Anderson;
    208:   snes->ops->reset          = SNESReset_NGMRES;

    210:   snes->usesnpc = PETSC_TRUE;
    211:   snes->usesksp = PETSC_FALSE;
    212:   snes->npcside = PC_RIGHT;

    214:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

    216:   PetscCall(PetscNew(&ngmres));
    217:   snes->data    = (void *)ngmres;
    218:   ngmres->msize = 30;

    220:   PetscCall(SNESParametersInitialize(snes));
    221:   PetscObjectParameterSetDefault(snes, max_funcs, 30000);
    222:   PetscObjectParameterSetDefault(snes, max_its, 10000);

    224:   PetscCall(SNESGetLineSearch(snes, &linesearch));
    225:   if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC));

    227:   ngmres->additive_linesearch = NULL;
    228:   ngmres->approxfunc          = PETSC_FALSE;
    229:   ngmres->restart_type        = SNES_NGMRES_RESTART_NONE;
    230:   ngmres->restart_it          = 2;
    231:   ngmres->restart_periodic    = 30;
    232:   ngmres->gammaA              = 2.0;
    233:   ngmres->gammaC              = 2.0;
    234:   ngmres->deltaB              = 0.9;
    235:   ngmres->epsilonB            = 0.1;

    237:   ngmres->andersonBeta = 1.0;
    238:   PetscFunctionReturn(PETSC_SUCCESS);
    239: }