Actual source code: viss.c

      1: #include <../src/snes/impls/vi/ss/vissimpl.h>

      3: /*@
      4:   SNESVIComputeMeritFunction - Evaluates the merit function for the mixed complementarity problem.

      6:   Input Parameter:
      7: . phi - the `Vec` holding the evaluation of the semismooth function

      9:   Output Parameters:
     10: + merit   - the merit function 1/2 ||phi||^2
     11: - phinorm - the two-norm of the vector, ||phi||

     13:   Level: developer

     15: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONSSLS`, `SNESVIComputeFunction()`
     16: @*/
     17: PetscErrorCode SNESVIComputeMeritFunction(Vec phi, PetscReal *merit, PetscReal *phinorm)
     18: {
     19:   PetscFunctionBegin;
     20:   PetscCall(VecNormBegin(phi, NORM_2, phinorm));
     21:   PetscCall(VecNormEnd(phi, NORM_2, phinorm));
     22:   *merit = 0.5 * (*phinorm) * (*phinorm);
     23:   PetscFunctionReturn(PETSC_SUCCESS);
     24: }

     26: static inline PetscScalar Phi(PetscScalar a, PetscScalar b)
     27: {
     28:   return a + b - PetscSqrtScalar(a * a + b * b);
     29: }

     31: static inline PetscScalar DPhi(PetscScalar a, PetscScalar b)
     32: {
     33:   if ((PetscAbsScalar(a) >= 1.e-6) || (PetscAbsScalar(b) >= 1.e-6)) return 1.0 - a / PetscSqrtScalar(a * a + b * b);
     34:   else return .5;
     35: }

     37: /*@
     38:   SNESVIComputeFunction - Provides the function that reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear
     39:   equations in semismooth form.

     41:   Input Parameters:
     42: + snes   - the `SNES` context
     43: . X      - current iterate
     44: - functx - user defined function context

     46:   Output Parameter:
     47: . phi - the evaluation of semismooth function at `X`

     49:   Level: developer

     51: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONSSLS`, `SNESVIComputeMeritFunction()`
     52: @*/
     53: PetscErrorCode SNESVIComputeFunction(SNES snes, Vec X, Vec phi, void *functx)
     54: {
     55:   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;
     56:   Vec                Xl = snes->xl, Xu = snes->xu, F = snes->vec_func;
     57:   PetscScalar       *phi_arr, *f_arr, *l, *u;
     58:   const PetscScalar *x_arr;
     59:   PetscInt           i, nlocal;

     61:   PetscFunctionBegin;
     62:   PetscCall((*vi->computeuserfunction)(snes, X, F, functx));
     63:   PetscCall(VecGetLocalSize(X, &nlocal));
     64:   PetscCall(VecGetArrayRead(X, &x_arr));
     65:   PetscCall(VecGetArray(F, &f_arr));
     66:   PetscCall(VecGetArray(Xl, &l));
     67:   PetscCall(VecGetArray(Xu, &u));
     68:   PetscCall(VecGetArray(phi, &phi_arr));

     70:   for (i = 0; i < nlocal; i++) {
     71:     if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
     72:       phi_arr[i] = f_arr[i];
     73:     } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
     74:       phi_arr[i] = -Phi(u[i] - x_arr[i], -f_arr[i]);
     75:     } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
     76:       phi_arr[i] = Phi(x_arr[i] - l[i], f_arr[i]);
     77:     } else if (l[i] == u[i]) {
     78:       phi_arr[i] = l[i] - x_arr[i];
     79:     } else { /* both bounds on variable */
     80:       phi_arr[i] = Phi(x_arr[i] - l[i], -Phi(u[i] - x_arr[i], -f_arr[i]));
     81:     }
     82:   }

     84:   PetscCall(VecRestoreArrayRead(X, &x_arr));
     85:   PetscCall(VecRestoreArray(F, &f_arr));
     86:   PetscCall(VecRestoreArray(Xl, &l));
     87:   PetscCall(VecRestoreArray(Xu, &u));
     88:   PetscCall(VecRestoreArray(phi, &phi_arr));
     89:   PetscFunctionReturn(PETSC_SUCCESS);
     90: }

     92: /*
     93:    SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the
     94:                                           the semismooth jacobian.
     95: */
     96: static PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes, Vec X, Vec F, Mat jac, Vec Da, Vec Db)
     97: {
     98:   PetscScalar *l, *u, *x, *f, *da, *db, da1, da2, db1, db2;
     99:   PetscInt     i, nlocal;

    101:   PetscFunctionBegin;
    102:   PetscCall(VecGetArray(X, &x));
    103:   PetscCall(VecGetArray(F, &f));
    104:   PetscCall(VecGetArray(snes->xl, &l));
    105:   PetscCall(VecGetArray(snes->xu, &u));
    106:   PetscCall(VecGetArray(Da, &da));
    107:   PetscCall(VecGetArray(Db, &db));
    108:   PetscCall(VecGetLocalSize(X, &nlocal));

    110:   for (i = 0; i < nlocal; i++) {
    111:     if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
    112:       da[i] = 0;
    113:       db[i] = 1;
    114:     } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
    115:       da[i] = DPhi(u[i] - x[i], -f[i]);
    116:       db[i] = DPhi(-f[i], u[i] - x[i]);
    117:     } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
    118:       da[i] = DPhi(x[i] - l[i], f[i]);
    119:       db[i] = DPhi(f[i], x[i] - l[i]);
    120:     } else if (l[i] == u[i]) { /* fixed variable */
    121:       da[i] = 1;
    122:       db[i] = 0;
    123:     } else { /* upper and lower bounds on variable */
    124:       da1   = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i]));
    125:       db1   = DPhi(-Phi(u[i] - x[i], -f[i]), x[i] - l[i]);
    126:       da2   = DPhi(u[i] - x[i], -f[i]);
    127:       db2   = DPhi(-f[i], u[i] - x[i]);
    128:       da[i] = da1 + db1 * da2;
    129:       db[i] = db1 * db2;
    130:     }
    131:   }

    133:   PetscCall(VecRestoreArray(X, &x));
    134:   PetscCall(VecRestoreArray(F, &f));
    135:   PetscCall(VecRestoreArray(snes->xl, &l));
    136:   PetscCall(VecRestoreArray(snes->xu, &u));
    137:   PetscCall(VecRestoreArray(Da, &da));
    138:   PetscCall(VecRestoreArray(Db, &db));
    139:   PetscFunctionReturn(PETSC_SUCCESS);
    140: }

    142: /*
    143:    SNESVIComputeJacobian - Computes the jacobian of the semismooth function.The Jacobian for the semismooth function is an element of the B-subdifferential of the Fischer-Burmeister function for complementarity problems.

    145:    Input Parameters:
    146: .  Da       - Diagonal shift vector for the semismooth jacobian.
    147: .  Db       - Row scaling vector for the semismooth jacobian.

    149:    Output Parameters:
    150: .  jac      - semismooth jacobian
    151: .  jac_pre  - optional preconditioning matrix

    153:    Note:
    154:    The semismooth jacobian matrix is given by
    155:    jac = Da + Db*jacfun
    156:    where Db is the row scaling matrix stored as a vector,
    157:          Da is the diagonal perturbation matrix stored as a vector
    158:    and   jacfun is the jacobian of the original nonlinear function.
    159: */
    160: static PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre, Vec Da, Vec Db)
    161: {
    162:   /* Do row scaling  and add diagonal perturbation */
    163:   PetscFunctionBegin;
    164:   PetscCall(MatDiagonalScale(jac, Db, NULL));
    165:   PetscCall(MatDiagonalSet(jac, Da, ADD_VALUES));
    166:   if (jac != jac_pre) { /* If jac and jac_pre are different */
    167:     PetscCall(MatDiagonalScale(jac_pre, Db, NULL));
    168:     PetscCall(MatDiagonalSet(jac_pre, Da, ADD_VALUES));
    169:   }
    170:   PetscFunctionReturn(PETSC_SUCCESS);
    171: }

    173: // PetscClangLinter pragma disable: -fdoc-sowing-chars
    174: /*
    175:   SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi.

    177:   Input Parameters:
    178:    phi - semismooth function.
    179:    H   - semismooth jacobian

    181:   Output Parameter:
    182:    dpsi - merit function gradient

    184:   Note:
    185:   The merit function gradient is computed as follows
    186:   dpsi = H^T*phi
    187: */
    188: static PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi)
    189: {
    190:   PetscFunctionBegin;
    191:   PetscCall(MatMultTranspose(H, phi, dpsi));
    192:   PetscFunctionReturn(PETSC_SUCCESS);
    193: }

    195: static PetscErrorCode SNESSolve_VINEWTONSSLS(SNES snes)
    196: {
    197:   SNES_VINEWTONSSLS   *vi = (SNES_VINEWTONSSLS *)snes->data;
    198:   PetscInt             maxits, i, lits;
    199:   SNESLineSearchReason lssucceed;
    200:   PetscReal            gnorm, xnorm = 0, ynorm;
    201:   Vec                  Y, X, F;
    202:   KSPConvergedReason   kspreason;
    203:   DM                   dm;
    204:   DMSNES               sdm;

    206:   PetscFunctionBegin;
    207:   PetscCall(SNESGetDM(snes, &dm));
    208:   PetscCall(DMGetDMSNES(dm, &sdm));

    210:   vi->computeuserfunction   = sdm->ops->computefunction;
    211:   sdm->ops->computefunction = SNESVIComputeFunction;

    213:   snes->numFailures            = 0;
    214:   snes->numLinearSolveFailures = 0;
    215:   snes->reason                 = SNES_CONVERGED_ITERATING;

    217:   maxits = snes->max_its;  /* maximum number of iterations */
    218:   X      = snes->vec_sol;  /* solution vector */
    219:   F      = snes->vec_func; /* residual vector */
    220:   Y      = snes->work[0];  /* work vectors */

    222:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    223:   snes->iter = 0;
    224:   snes->norm = 0.0;
    225:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

    227:   PetscCall(SNESVIProjectOntoBounds(snes, X));
    228:   PetscCall(SNESComputeFunction(snes, X, vi->phi));
    229:   if (snes->domainerror) {
    230:     snes->reason              = SNES_DIVERGED_FUNCTION_DOMAIN;
    231:     sdm->ops->computefunction = vi->computeuserfunction;
    232:     PetscFunctionReturn(PETSC_SUCCESS);
    233:   }
    234:   /* Compute Merit function */
    235:   PetscCall(SNESVIComputeMeritFunction(vi->phi, &vi->merit, &vi->phinorm));

    237:   PetscCall(VecNormBegin(X, NORM_2, &xnorm)); /* xnorm <- ||x||  */
    238:   PetscCall(VecNormEnd(X, NORM_2, &xnorm));
    239:   SNESCheckFunctionNorm(snes, vi->merit);

    241:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    242:   snes->norm = vi->phinorm;
    243:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
    244:   PetscCall(SNESLogConvergenceHistory(snes, vi->phinorm, 0));

    246:   /* test convergence */
    247:   PetscCall(SNESConverged(snes, 0, 0.0, 0.0, vi->phinorm));
    248:   PetscCall(SNESMonitor(snes, 0, vi->phinorm));
    249:   if (snes->reason) {
    250:     sdm->ops->computefunction = vi->computeuserfunction;
    251:     PetscFunctionReturn(PETSC_SUCCESS);
    252:   }

    254:   for (i = 0; i < maxits; i++) {
    255:     /* Call general purpose update function */
    256:     PetscTryTypeMethod(snes, update, snes->iter);

    258:     /* Solve J Y = Phi, where J is the semismooth jacobian */

    260:     /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
    261:     sdm->ops->computefunction = vi->computeuserfunction;
    262:     PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
    263:     SNESCheckJacobianDomainerror(snes);
    264:     sdm->ops->computefunction = SNESVIComputeFunction;

    266:     /* Get the diagonal shift and row scaling vectors */
    267:     PetscCall(SNESVIComputeBsubdifferentialVectors(snes, X, F, snes->jacobian, vi->Da, vi->Db));
    268:     /* Compute the semismooth jacobian */
    269:     PetscCall(SNESVIComputeJacobian(snes->jacobian, snes->jacobian_pre, vi->Da, vi->Db));
    270:     /* Compute the merit function gradient */
    271:     PetscCall(SNESVIComputeMeritFunctionGradient(snes->jacobian, vi->phi, vi->dpsi));
    272:     PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
    273:     PetscCall(KSPSolve(snes->ksp, vi->phi, Y));
    274:     PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason));

    276:     if (kspreason < 0) {
    277:       if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
    278:         PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n", snes->iter, snes->numLinearSolveFailures));
    279:         snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
    280:         break;
    281:       }
    282:     }
    283:     PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
    284:     snes->linear_its += lits;
    285:     PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
    286:     /*
    287:     if (snes->ops->precheck) {
    288:       PetscBool changed_y = PETSC_FALSE;
    289:       PetscUseTypeMethod(snes,precheck ,X,Y,snes->precheck,&changed_y);
    290:     }

    292:     if (PetscLogPrintInfo) PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W));
    293:     */
    294:     /* Compute a (scaled) negative update in the line search routine:
    295:          Y <- X - lambda*Y
    296:        and evaluate G = function(Y) (depends on the line search).
    297:     */
    298:     PetscCall(VecCopy(Y, snes->vec_sol_update));
    299:     ynorm = 1;
    300:     gnorm = vi->phinorm;
    301:     PetscCall(SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y));
    302:     PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed));
    303:     PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm));
    304:     PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)vi->phinorm, (double)gnorm, (double)ynorm, (int)lssucceed));
    305:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
    306:     if (snes->domainerror) {
    307:       snes->reason              = SNES_DIVERGED_FUNCTION_DOMAIN;
    308:       sdm->ops->computefunction = vi->computeuserfunction;
    309:       PetscFunctionReturn(PETSC_SUCCESS);
    310:     }
    311:     if (lssucceed) {
    312:       if (++snes->numFailures >= snes->maxFailures) {
    313:         PetscBool ismin;
    314:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
    315:         PetscCall(SNESVICheckLocalMin_Private(snes, snes->jacobian, vi->phi, X, gnorm, &ismin));
    316:         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
    317:         break;
    318:       }
    319:     }
    320:     /* Update function and solution vectors */
    321:     vi->phinorm = gnorm;
    322:     vi->merit   = 0.5 * vi->phinorm * vi->phinorm;
    323:     /* Monitor convergence */
    324:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    325:     snes->iter  = i + 1;
    326:     snes->norm  = vi->phinorm;
    327:     snes->xnorm = xnorm;
    328:     snes->ynorm = ynorm;
    329:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
    330:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
    331:     /* Test for convergence, xnorm = || X || */
    332:     if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm));
    333:     PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, vi->phinorm));
    334:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
    335:     if (snes->reason) break;
    336:   }
    337:   sdm->ops->computefunction = vi->computeuserfunction;
    338:   PetscFunctionReturn(PETSC_SUCCESS);
    339: }

    341: static PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
    342: {
    343:   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;

    345:   PetscFunctionBegin;
    346:   PetscCall(SNESSetUp_VI(snes));
    347:   PetscCall(VecDuplicate(snes->vec_sol, &vi->dpsi));
    348:   PetscCall(VecDuplicate(snes->vec_sol, &vi->phi));
    349:   PetscCall(VecDuplicate(snes->vec_sol, &vi->Da));
    350:   PetscCall(VecDuplicate(snes->vec_sol, &vi->Db));
    351:   PetscCall(VecDuplicate(snes->vec_sol, &vi->z));
    352:   PetscCall(VecDuplicate(snes->vec_sol, &vi->t));
    353:   PetscFunctionReturn(PETSC_SUCCESS);
    354: }

    356: static PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
    357: {
    358:   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;

    360:   PetscFunctionBegin;
    361:   PetscCall(SNESReset_VI(snes));
    362:   PetscCall(VecDestroy(&vi->dpsi));
    363:   PetscCall(VecDestroy(&vi->phi));
    364:   PetscCall(VecDestroy(&vi->Da));
    365:   PetscCall(VecDestroy(&vi->Db));
    366:   PetscCall(VecDestroy(&vi->z));
    367:   PetscCall(VecDestroy(&vi->t));
    368:   PetscFunctionReturn(PETSC_SUCCESS);
    369: }

    371: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(SNES snes, PetscOptionItems *PetscOptionsObject)
    372: {
    373:   PetscFunctionBegin;
    374:   PetscCall(SNESSetFromOptions_VI(snes, PetscOptionsObject));
    375:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES semismooth method options");
    376:   PetscOptionsHeadEnd();
    377:   PetscFunctionReturn(PETSC_SUCCESS);
    378: }

    380: /*MC
    381:       SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method

    383:    Options Database Keys:
    384: +   -snes_type <vinewtonssls,vinewtonrsls> a semi-smooth solver, a reduced space active set method
    385: -   -snes_vi_monitor - prints the number of active constraints at each iteration.

    387:    Level: beginner

    389:    Notes:
    390:    This family of algorithms is much like an interior point method.

    392:    The reduced space active set solvers `SNESVINEWTONRSLS` provide an alternative approach that does not result in extremely ill-conditioned linear systems

    394:    See {cite}`munson.facchinei.ea:semismooth` and {cite}`benson2006flexible`

    396: .seealso: [](ch_snes), `SNESVINEWTONRSLS`, `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONRSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`
    397: M*/
    398: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
    399: {
    400:   SNES_VINEWTONSSLS *vi;
    401:   SNESLineSearch     linesearch;

    403:   PetscFunctionBegin;
    404:   snes->ops->reset          = SNESReset_VINEWTONSSLS;
    405:   snes->ops->setup          = SNESSetUp_VINEWTONSSLS;
    406:   snes->ops->solve          = SNESSolve_VINEWTONSSLS;
    407:   snes->ops->destroy        = SNESDestroy_VI;
    408:   snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
    409:   snes->ops->view           = NULL;

    411:   snes->usesksp = PETSC_TRUE;
    412:   snes->usesnpc = PETSC_FALSE;

    414:   PetscCall(SNESGetLineSearch(snes, &linesearch));
    415:   if (!((PetscObject)linesearch)->type_name) {
    416:     PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
    417:     PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0));
    418:   }

    420:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

    422:   PetscCall(SNESParametersInitialize(snes));

    424:   PetscCall(PetscNew(&vi));
    425:   snes->data = (void *)vi;

    427:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI));
    428:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI));
    429:   PetscFunctionReturn(PETSC_SUCCESS);
    430: }